gavo.utils.pgsphere

1 """ 2 Bindings for the pgsphere libarary and psycopg2. 3 4 Basically, once per program run, you need to call preparePgSphere(connection), 5 and you're done. 6 7 All native representation is in rad. 8 """ 9 10 #c Copyright 2008-2019, the GAVO project 11 #c 12 #c This program is free software, covered by the GNU GPL. See the 13 #c COPYING file in the source distribution. 14 15 16 import math 17 import re 18 import tempfile 19 from cStringIO import StringIO 20 21 from gavo.utils import codetricks 22 from gavo.utils import excs 23 from gavo.utils import mathtricks 24 from gavo.utils import misctricks 25 from gavo.utils import texttricks 26 from gavo.utils.fitstools import pyfits 27 from gavo.utils.mathtricks import DEG 28 29 30 _TRAILING_ZEROES = re.compile("0+(\s|$)")

31 -def removeTrailingZeroes(s):

32 """remove zeroes in front of whitespace or the string end. 33 34 This is used for cosmetics in STC-S strings. 35 36 >>> removeTrailingZeroes("23.3420 21.2 12.00000") 37 '23.342 21.2 12.' 38 """ 39 return _TRAILING_ZEROES.sub(r"\1", s)

40

41 42 -class TwoSBoxes(excs.ExecutiveAction):

43 """is raised when an SBox is constructed from center and size such that 44 it overlaps the pole. 45 """

46 - def __init__(self, box1, box2):

47 self.box1, self.box2 = box1, box2

48

49 50 -def _query(conn, query, pars=None):

51 c = conn.cursor() 52 c.execute(query, pars) 53 res = list(c) 54 return res

55

56 57 -class PgSAdapter(object):

58 """A base class for objects adapting pgSphere objects. 59 60 The all need a pgType attribute and two static methods 61 _adaptToPgSphere(obj) and _castFromPgSphere(value, cursor). 62 63 You must also define a sequence checkedAttributes; all attributes 64 metioned there must be equal for two adapted values to be equal (equality 65 testing here really is mainly for unit tests with hand-crafted values). 66 67 Also, all subclasses you should provide an asPoly returning a spherical 68 polygon. This is used when uploading VOTables with REGION columns. 69 70 In practice, you also want: 71 72 fromDALI(cls, daliSeq) -- a classmethod turning the DALI votable 73 representation (typically, float arrays in degrees) into a corresponding 74 object. 75 asDALI(self) -- the inverse of fromDALI 76 asSODA(self) -- returns a representation of the geometry as 77 used in SODA parameters. 78 """ 79 pgType = None 80

81 - def __eq__(self, other):

82 if self.__class__!=other.__class__: 83 return False 84 for attName in self.checkedAttributes: 85 if getattr(self, attName)!=getattr(other, attName): 86 return False 87 return True

88

89 - def __ne__(self, other):

90 return not self==other

91

92 - def __repr__(self):

93 return "<pgsphere %s>"%self.asSTCS("Unknown")

94

95 - def asPoly(self):

96 raise ValueError("%s objects cannot be turned into polygons."% 97 self.__class__)

98

99 - def asSODA(self):

100 """returns the "SODA-form" for this circle. 101 102 This is a string containing blank-separated float literals; the 103 floats are what asDALI returns, and hence this is built on top of 104 asDALI. We don't worry about the fact that the coordinates just 105 *might* be non-ICRS. 106 """ 107 return removeTrailingZeroes( 108 " ".join("%.10f"%v for v in self.asDALI()))

109

110 111 -class SPoint(PgSAdapter):

112 """A point on a sphere from pgSphere. 113 114 >>> SPoint(1, 1).asSTCS("ICRS") 115 'Position ICRS 57.2957795131 57.2957795131' 116 >>> SPoint.fromDegrees(1, -1).asSTCS("ICRS") 117 'Position ICRS 1. -1.' 118 """ 119 pgType = "spoint" 120 checkedAttributes = ["x", "y"] 121 pattern = re.compile(r"\s*$\s*([0-9.e-]+)\s*,\s*([0-9.e-]+)\s*$") 122

123 - def __init__(self, x, y):

124 self.x, self.y = float(x), float(y)

125

126 - def __repr__(self):

127 return "SPoint(%r, %r)"%(self.x, self.y)

128 129 @staticmethod

130 - def _adaptToPgSphere(spoint):

131 return AsIs("spoint '(%.10f,%.10f)'"%(spoint.x, spoint.y))

132 133 @classmethod

134 - def _castFromPgSphere(cls, value, cursor):

135 if value is not None: 136 return cls(*map(float, cls.pattern.match(value).groups()))

137 138 @classmethod

139 - def fromDegrees(cls, x, y):

140 if x is None or y is None: 141 return None 142 return cls(x*DEG, y*DEG)

143

144 - def asCooPair(self):

145 """returns this point as (long, lat) in degrees. 146 """ 147 return (self.x/DEG, self.y/DEG)

148

149 - def asSTCS(self, systemString):

150 return removeTrailingZeroes( 151 "Position %s %.10f %.10f"%(systemString, self.x/DEG, self.y/DEG))

152

153 - def asPgSphere(self):

154 return "spoint '(%.10f,%.10f)'"%(self.x, self.y)

155

156 - def p(self): # helps below

157 return "(%r, %r)"%(self.x, self.y)

158

159 - def asUnitSphereVec(self):

160 """returns self as a triple of cx, cy, cz on the unit sphere. 161 """ 162 return mathtricks.spherToCart(self.x, self.y)

163

164 - def asDALI(self):

165 return self.asCooPair()

166 167 @classmethod

168 - def fromDALI(cls, coos):

169 return cls.fromDegrees(*coos)

170

171 172 -class SCircle(PgSAdapter):

173 """A spherical circle from pgSphere. 174 175 The constructor accepts an SPoint center and a radius in rad. 176 """ 177 pgType = "scircle" 178 checkedAttributes = ["center", "radius"] 179 pattern = re.compile("<($[^)]*$)\s*,\s*([0-9.e-]+)>") 180

181 - def __init__(self, center, radius):

182 self.center, self.radius = center, float(radius)

183 184 @staticmethod

185 - def _adaptToPgSphere(sc):

186 return AsIs("scircle '< %s, %r >'"%(sc.center.p(), sc.radius))

187 188 @classmethod

189 - def _castFromPgSphere(cls, value, cursor):

190 if value is not None: 191 pt, radius = cls.pattern.match(value).groups() 192 return cls(SPoint._castFromPgSphere(pt, cursor), radius)

193

194 - def asDALI(self):

195 """returns the "DALI-form" for this circle. 196 197 This is an array containing the center coordinates and the radius 198 in degrees. 199 """ 200 return (self.center.x/DEG, self.center.y/DEG, self.radius/DEG)

201 202 @classmethod

203 - def fromDALI(cls, daliSeq):

204 """returns a circle from its DALI float sequence. 205 206 Any None in daliSeq makes this a None to help with VOTable null value 207 handling. 208 """ 209 if None in daliSeq: 210 return None 211 ra, dec, radius = [float(s) for s in daliSeq] 212 return cls(SPoint.fromDegrees(ra, dec), radius*DEG)

213

214 - def asSTCS(self, systemString):

215 return removeTrailingZeroes("Circle %s %s"%( 216 systemString, self.asSODA()))

217

218 - def asPgSphere(self):

219 return "scircle '< (%.10f, %.10f), %.10f >'"%( 220 self.center.x, self.center.y, self.radius)

221

222 - def asPoly(self, nSegments=32):

223 # approximate the circle with 32 line segments and don't worry about 224 # circles with radii larger than 90 degrees. 225 # We compute the circle around the north pole and then rotate 226 # the resulting points such that the center ends up at the 227 # circle's center. 228 r = math.sin(self.radius) 229 innerOffset = math.cos(self.radius) 230 rotationMatrix = mathtricks.getRotZ(math.pi/2-self.center.x).matMul( 231 mathtricks.getRotX(math.pi/2-self.center.y)) 232 233 points = [] 234 for i in range(nSegments): 235 angle = 2.*i*math.pi/nSegments 236 dx, dy = r*math.sin(angle), r*math.cos(angle) 237 points.append(SPoint( 238 *mathtricks.cartToSpher(rotationMatrix.vecMul((dx, dy, innerOffset))))) 239 return SPoly(points)

240

241 242 -class SPoly(PgSAdapter):

243 """A spherical polygon from pgSphere. 244 245 The constructor accepts a list points of SPoints. 246 """ 247 pgType = "spoly" 248 checkedAttributes = ["points"] 249 pattern = re.compile("$[^)]+$") 250

251 - def __init__(self, points):

252 self.points = tuple(points)

253 254 @staticmethod

255 - def _adaptToPgSphere(spoly):

256 return AsIs("spoly '{%s}'"%(", ".join(p.p() for p in spoly.points)))

257 258 @classmethod

259 - def _castFromPgSphere(cls, value, cursor):

260 if value is not None: 261 return cls([SPoint._castFromPgSphere(ptLit, cursor) 262 for ptLit in cls.pattern.findall(value)])

263

264 - def asDALI(self):

265 """returns the DALI form of this polygon. 266 267 This is a sequence of floats of the vertex coordinates in degrees. 268 """ 269 res = [] 270 for p in self.points: 271 res.extend([p.x/DEG, p.y/DEG]) 272 return res

273 274 @classmethod

275 - def fromDALI(cls, daliSeq):

276 """returns a polygon from a DALI float-sequence 277 278 This returns None for empty daliSeqs to help with VOTable NULL value 279 handling. 280 """ 281 if not daliSeq: 282 return None 283 return cls([SPoint.fromDegrees(*tuple(p)) 284 for p in misctricks.grouped(2, daliSeq)])

285

286 - def asCooPairs(self):

287 """returns the vertices as a sequence of (long, lat) pairs in 288 degrees. 289 290 This form is required by some functions from base.coords. 291 """ 292 return [p.asCooPair() for p in self.points]

293

294 - def asSTCS(self, systemString):

295 return removeTrailingZeroes("Polygon %s %s"%(systemString, 296 self.asSODA()))

297

298 - def asPgSphere(self):

299 return "spoly '{%s}'"%(",".join("(%.10f,%.10f)"%(p.x, p.y) 300 for p in self.points))

301

302 - def asPoly(self):

303 return self

304

305 - def asSMoc(self, order=6, inclusive=True):

306 """returns an SMoc instance for this polygon. 307 308 If inclusive is False, do not include healpixes that lie 309 310 order is the maximum order of the moc returned. 311 """ 312 import healpy 313 pixels = healpy.query_polygon(vertices=[ 314 p.asUnitSphereVec() for p in self.points], 315 nside=2**order, nest=True) 316 return SMoc.fromCells(order, pixels)

317

318 319 -class SBox(PgSAdapter):

320 """A spherical box from pgSphere. 321 322 The constructor accepts the two corner points. 323 """ 324 pgType = "sbox" 325 checkedAttributes = ["corner1", "corner2"] 326 pattern = re.compile("$[^()]+$") 327

328 - def __init__(self, corner1, corner2):

329 self.corner1, self.corner2 = corner1, corner2

330 331 @staticmethod

332 - def _adaptToPgSphere(sbox):

333 return AsIs("sbox '(%s, %s)'"%(sbox.corner1.p(), sbox.corner2.p()))

334 335 @classmethod

336 - def _castFromPgSphere(cls, value, cursor):

337 if value is not None: 338 return cls(*[SPoint._castFromPgSphere(ptLit, cursor) 339 for ptLit in cls.pattern.findall(value)])

340 341 @classmethod

342 - def fromSIAPPars(cls, ra, dec, raSize, decSize):

343 """returns an SBox corresponding to what SIAP passes in. 344 345 In particular, all values are in degrees, and a cartesian projection 346 is assumed. 347 348 This is for use with SIAP and tries to partially implement that silly 349 prescription of "folding" over at the poles. If that happens, 350 a TwoSBoxes exception is raised. It contains two SBoxes that 351 should be ORed. I agree that sucks. Let's fix SIAP. 352 """ 353 if 90-abs(dec)<0.1: # Special handling at the pole 354 raSize = 360 355 else: 356 raSize = raSize/math.cos(dec*DEG) 357 decSize = abs(decSize) # inhibit auto swapping of points 358 minRA, maxRA = ra-raSize/2., ra+raSize/2. 359 bottom, top = dec-decSize/2., dec+decSize/2. 360 # folding over at the poles: raise an exception with two boxes, 361 # and let upstream handle it. Foldover on both poles is not supported. 362 # All this isn't really thought out and probably doesn't work in 363 # many interesting cases. 364 # I hate that folding over. 365 if bottom<-90 and top>90: 366 raise ValueError("Cannot fold over at both poles") 367 elif bottom<-90: 368 raise TwoSBoxes( 369 cls( 370 SPoint.fromDegrees(minRA, -90), 371 SPoint.fromDegrees(maxRA, top)), 372 cls( 373 SPoint.fromDegrees(180+minRA, -90), 374 SPoint.fromDegrees(180+maxRA, top))) 375 elif top>90: 376 raise TwoSBoxes( 377 cls( 378 SPoint.fromDegrees(minRA, bottom), 379 SPoint.fromDegrees(maxRA, 90)), 380 cls( 381 SPoint.fromDegrees(180+minRA, bottom), 382 SPoint.fromDegrees(180+maxRA, 90))) 383 return cls(SPoint.fromDegrees(minRA, bottom), 384 SPoint.fromDegrees(maxRA, top))

385

386 - def asSTCS(self, systemString):

387 return removeTrailingZeroes("PositionInterval %s %s %s"%(systemString, 388 "%.10f %.10f"%(self.corner1.x/DEG, self.corner1.y/DEG), 389 "%.10f %.10f"%(self.corner2.x/DEG, self.corner2.y/DEG)))

390

391 - def asPoly(self):

392 x1, y1 = self.corner1.x, self.corner1.y 393 x2, y2 = self.corner2.x, self.corner2.y 394 minX, maxX = min(x1, x2), max(x1, x2) 395 minY, maxY = min(y1, y2), max(y1, y2) 396 return SPoly(( 397 SPoint(minX, minY), 398 SPoint(minX, maxY), 399 SPoint(maxX, maxY), 400 SPoint(maxX, minY)))

401 402 @classmethod

403 - def fromDALI(cls, daliSeq):

404 # see asDALI() 405 return cls( 406 SPoint.fromDALI(daliSeq[:2]), 407 SPoint.fromDALI(daliSeq[2:]))

408

409 - def asDALI(self):

410 # we're cheating a bit here; there's no official DALI representation 411 # yet. This is for our private xtype x:box 412 return self.corner1.asDALI()+self.corner2.asDALI()

413 414 415 # temporary measure to shut up astropy's configuration parser. 416 import __builtin__ 417 __builtin__._ASTROPY_SETUP_ = True

418 419 420 -class SMoc(PgSAdapter):

421 """a MOC with a pgsphere interface. 422 423 The default constructor accepts a pymoc.MOC instance, which is 424 accessible as the moc attribute. The database interface uses 425 the ASCII serialisation, for which there's the fromASCII constructor. 426 """ 427 pgType = "smoc" 428 checkedAttributes = ["moc"] 429

430 - def __init__(self, moc):

431 self.moc = moc

432 433 _orderHeadRE = re.compile(r"(\d+)/") 434 _rangeRE = re.compile(r"\s*(\d+)(?:-(\d+))?\s*$") 435 436 @classmethod

437 - def _parseCells(cls, cellLiteral, firstPos):

438 """returns a sequence of cells from a MOC cell literal. 439 440 firstPos is the string position of the beginning of cellLiteral 441 for error messages. 442 """ 443 curPos = 0 444 cells = [] 445 446 for item in cellLiteral.split(","): 447 mat = cls._rangeRE.match(item) 448 if not mat: 449 raise ValueError("MOC literal syntax error at char %s"% 450 (firstPos+curPos)) 451 452 if mat.group(2) is None: 453 cells.append(int(mat.group(1))) 454 else: 455 cells.extend(range(int(mat.group(1)), int(mat.group(2))+1)) 456 457 curPos += len(item)+1 458 459 return cells

460 461 @classmethod

462 - def fromASCII(cls, literal):

463 """returns an SMoc from a quasi-standard ASCII serialisation. 464 """ 465 import pymoc 466 467 # we do the parsing ourselves -- the pymoc interface is too clumsy, 468 # and the rigidity of the parser doesn't look good. 469 seps = list(cls._orderHeadRE.finditer(literal))+[re.search("$", literal)] 470 if len(seps)==1: 471 raise ValueError("No order separator visible in MOC literal '%s'"% 472 texttricks.makeEllipsis(literal, 40)) 473 if not re.match(r"\s*$", literal[:seps[0].start()]): 474 raise ValueError("MOC literal '%s' does not start with order spec"% 475 texttricks.makeEllipsis(literal, 40)) 476 477 moc = pymoc.MOC() 478 479 for openMat, closeMat in codetricks.iterRanges(seps): 480 order = int(openMat.group(1)) 481 cells = cls._parseCells(literal[openMat.end():closeMat.start()], 482 openMat.end()) 483 moc.add(order, cells) 484 485 return cls(moc)

486 487 @classmethod

488 - def fromCells(cls, order, pixels):

489 """returns a SMoc instance from a collection of pixels at order. 490 """ 491 import pymoc 492 493 moc = pymoc.MOC(order=order, cells=pixels) 494 moc.normalize() 495 return cls(moc)

496 497 @staticmethod

498 - def _formatASCIIRange(minCell, maxCell):

499 """returns a cell literal for a MOC. 500 """ 501 if minCell==maxCell: 502 return str(minCell) 503 else: 504 return "%d-%d"%(minCell, maxCell)

505

506 - def asASCII(self):

507 """returns an ascii serialisation of this MOC. 508 """ 509 # this is essentially the pymoc code, but again saving the file 510 # interface that we don't want here. 511 parts = [] 512 for order, cells in self.moc: 513 ranges = [] 514 rmin = rmax = None 515 516 for cell in sorted(cells): 517 if rmin is None: 518 rmin = rmax = cell 519 elif rmax==cell-1: 520 rmax = cell 521 else: 522 ranges.append(self._formatASCIIRange(rmin, rmax)) 523 rmin = rmax = cell 524 525 ranges.append(self._formatASCIIRange(rmin, rmax)) 526 parts.append("%d/%s"%(order, ",".join(ranges))) 527 528 return " ".join(parts)

529 530 @classmethod

531 - def fromFITS(cls, literal):

532 """returns an SMoc from a string containing a FITS-serialised MOC. 533 """ 534 from pymoc.io.fits import read_moc_fits_hdu 535 import pymoc 536 537 moc = pymoc.MOC() 538 read_moc_fits_hdu(moc, 539 pyfits.open(StringIO(literal))[1]) 540 return cls(moc)

541 542 @staticmethod

543 - def _adaptToPgSphere(smoc):

544 return AsIs("smoc '%s'"%(smoc.asASCII()))

545 546 @classmethod

547 - def _castFromPgSphere(cls, value, cursor):

548 if value is not None: 549 return cls.fromASCII(value)

550

551 - def asPoly(self):

552 raise TypeError("MOCs cannot be represented as polygons")

553

554 - def asSTCS(self, frame):

555 # no STCS for MOCs, but this is really just for old-style VOTable 556 # serialisation, so let's cheat 557 return "MOC "+self.asASCII()

558

559 - def asSMoc(self, order=6):

560 """returns a copy of self, normalised for order. 561 """ 562 moc = self.moc.copy() 563 moc.normalize(order) 564 return self.__class__(moc)

565

566 - def getPlot(self, **kwargs):

567 """returns a png string with a plot visualising this moc. 568 """ 569 from pymoc.util.plot import plot_moc 570 with tempfile.NamedTemporaryFile(suffix=".png") as f: 571 plot_moc(self.moc, filename=f.name, projection="moll", **kwargs) 572 return f.read()

573

574 - def asFITS(self):

575 """returns a standard FITS (table) representation of this MOC. 576 """ 577 from pymoc.io.fits import write_moc_fits 578 579 with tempfile.NamedTemporaryFile(suffix=".fits") as f: 580 write_moc_fits(self.moc, f) 581 f.seek(0) 582 return f.read()

583

584 - def asDALI(self):

585 """returns the string representation of this MOC. 586 587 This isn't what DALI actually says at this point, but we suspect 588 it will say that at some point. 589 """ 590 return self.asASCII()

591 592 @classmethod

593 - def fromDALI(cls, literal):

594 """returns an SMoc from a MOC string. 595 596 (see asDALI) 597 """ 598 return cls.fromASCII(literal)

599

600 - def asNUNIQs(self):

601 """returns a list of integers usable as nuniqs. 602 """ 603 res = [] 604 for order, cells in self.moc: 605 shiftedOrder = 4 * (4 ** order) 606 res.extend(pix+shiftedOrder for pix in cells) 607 return res

608 609 610 try: 611 import psycopg2 612 from psycopg2.extensions import (register_adapter, AsIs, register_type, 613 new_type) 614 615 _getPgSClass = codetricks.buildClassResolver(PgSAdapter, globals().values(), 616 key=lambda obj: obj.pgType, default=PgSAdapter)

617 618 619 - def preparePgSphere(conn):

620 if hasattr(psycopg2, "_pgsphereLoaded"): 621 return 622 try: 623 oidmap = _query(conn, 624 "SELECT typname, oid" 625 " FROM pg_type" 626 " WHERE typname ~ " 627 " '^s(point|trans|circle|line|ellipse|poly|path|box|moc)'") 628 for typeName, oid in oidmap: 629 cls = _getPgSClass(typeName) 630 if cls is not PgSAdapter: # base class is null value 631 register_adapter(cls, cls._adaptToPgSphere) 632 register_type( 633 new_type((oid,), "spoint", cls._castFromPgSphere)) 634 psycopg2._pgsphereLoaded = True 635 conn.commit() 636 except: 637 psycopg2._pgsphereLoaded = False

638 639 except ImportError: 640 raise 641 # psycopg2 not installed. Since preparsePgSphere can only be 642 # called from code depending on psycopg2, there's no harm if 643 # we don't define it. 644 pass

645 646 647 -def _test():

648 import doctest, pgsphere 649 doctest.testmod(pgsphere)

650 651 652 if __name__=="__main__": 653 _test() 654

Source Code for Module gavo.utils.pgsphere