gavo.stc.tapstc

1 """ 2 The subset of STC proposed by the TAP spec. 3 4 Use this rather than the full-blown STC library for TAP and friends. 5 TAP's semi-sanitation of STC needs some special handling anyway, 6 and this is much faster than the full-blown mess. 7 """ 8 9 #c Copyright 2008-2019, the GAVO project 10 #c 11 #c This program is free software, covered by the GNU GPL. See the 12 #c COPYING file in the source distribution. 13 14 15 from __future__ import print_function 16 17 from gavo import utils 18 from gavo.stc import common 19 from gavo.stc import stcsast 20 from gavo.utils import pgsphere 21 22 23 TAP_SYSTEMS = set( 24 ['ICRS', 'FK4', 'FK5', 'GALACTIC', 'RELOCATABLE', 'UNKNOWN', '', "BROKEN", 25 "UNKNOWNFrame"]) 26 27 # universally ignored 28 TAP_REFPOS = set( 29 ['BARYCENTER', 'GEOCENTER', 'HELIOCENTER', 'LSR', 'TOPOCENTER', 30 'RELOCATABLE', 'UNKNOWNREFPOS']) 31 32 # only SPHERICAL2 supported, all others raise errors 33 TAP_FLAVORS = set( 34 ["CARTESIAN2", "CARTESIAN3", "SPHERICAL2"]) 35 36 37 ################## transformations between TAP STC reference systems 38 39 UNIVERSALLY_COMPATIBLE = set(['RELOCATABLE', 'UNKNOWN', '', "BROKEN", None]) 40 41 TRANSFORMS = { 42 # From "ICRS" (really, FK5 J2000, and the lo-prec) to... 43 'FK4': (1.5651864333666516, -0.0048590552804904244, -1.5763681043529187), 44 'FK5': None, 45 'ICRS': None, 46 'GALACTIC': (1.3463560974407338, -1.0973190018372752, 0.57477052472873258), 47 }

48 49 50 -def _getEulersFor(frame):

51 if not frame in TRANSFORMS: 52 raise common.STCValueError("Unknown reference frame: %s"%frame) 53 return TRANSFORMS[frame]

54

55 56 -def getPGSphereTrafo(fromSys, toSys):

57 """returns a pgsphere expression fragment to turn a pgsphere geometry 58 from fromSys to toSys. 59 60 fromSys and toSys are system designations like for TAP. 61 62 If the systems are "compatible" (i.e., no transformation is deemed necessary), 63 None is returned. An STCValueError is raised for incomprehensible system 64 specifications. 65 """ 66 if (fromSys in UNIVERSALLY_COMPATIBLE 67 or toSys in UNIVERSALLY_COMPATIBLE): 68 return None 69 if fromSys=='ICRS': 70 template = "+strans(%.12f,%.12f,%.12f)" 71 angles = _getEulersFor(toSys) 72 elif toSys=='ICRS': 73 angles = _getEulersFor(fromSys) 74 template = "-strans(%.12f,%.12f,%.12f)" 75 else: 76 t1 = getPGSphereTrafo(fromSys, 'ICRS') 77 t2 = getPGSphereTrafo('ICRS', toSys) 78 return "%s%s"%(t1 or "", t2 or "") 79 if angles: 80 return template%angles 81 return None

82

83 84 -def getTAPSTC(stcInstance):

85 """returns a tap system identifier for an STC AstroSystem. 86 87 This is stc.astroSystem.spaceFrame.refFrame if existing and 88 TAP-defined, UNKNOWN otherwise. 89 """ 90 rf = None 91 if stcInstance.astroSystem and stcInstance.astroSystem.spaceFrame: 92 rf = stcInstance.astroSystem.spaceFrame.refFrame 93 if rf not in TAP_SYSTEMS: 94 return "UNKNOWN" 95 return rf

96

97 98 @utils.memoized 99 -def getSTCForTAP(tapIdentifier):

100 """returns an stc AST for a tap reference system identifier. 101 102 The tapIdentifier is any string in which the first item is the reference 103 system. Everything else is ignored (this is because it seems someone 104 intended additional specs like "BARYCENTER" to be legal, although 105 there really is nothing we can do about them). 106 """ 107 if tapIdentifier: 108 tapIdentifier = utils.identifierPattern.findall(tapIdentifier)[0] 109 if tapIdentifier in ["BROKEN", '', "UNKNOWN"]: 110 tapIdentifier = "UNKNOWNFrame" 111 ast = stcsast.parseSTCS("Position %s"%tapIdentifier) 112 if tapIdentifier=='BROKEN': 113 ast.broken = True 114 return ast

115

116 117 118 ############################# TAP simplified STC-S 119 # The simplified STC-S is a simple regular grammar for the position specs 120 # plus a context-free part for set operations. The regular part we 121 # do with REs, for the rest there's a simple recursive descent parser. 122 # 123 # From the literal we either create a pgsphere geometry ready for ingestion 124 # or, when operators enter, an GeomExpr object. Since at least pgsphere 125 # does not support geometric operators, this must be handled in the morph code. 126 # To make this clear, its flatten method just raises an Exception. 127 # 128 # The regions expressible in pgsphere are returned as pgsphre objects. 129 # This is because I feel all this should only be used for ingesting data 130 # ("table upload") and thus carrying around the frame is pointless; 131 # you can, however, retrieve a guess for the frame from the returned 132 # object's cooSys attribute (which is convenient for the morphers). 133 134 135 -class GeomExpr(object):

136 """a sentinel object to be processed by morphers. 137 """

138 - def __init__(self, frame, operator, operands):

139 self.cooSys = frame 140 self.operator, self.operands = operator.upper(), operands

141

142 - def flatten(self):

143 raise common.STCSParseError("Cannot serialize STC-S. Did you use" 144 " Union or Intersection outside of CONTAINS or INTERSECTS?")

145

146 - def _flatLogic(self, template, operand):

147 if isinstance(operand, GeomExpr): 148 return operand.asLogic(template) 149 else: 150 return template%operand.asPgSphere()

151

152 - def asLogic(self, template):

153 if self.operator=="UNION": 154 logOp = " OR " 155 elif self.operator=="INTERSECTION": 156 logOp = " AND " 157 elif self.operator=="NOT": 158 return "NOT (%s)"%self._flatLogic(template, self.operands[0]) 159 else: 160 raise NotImplementedError("No logic for operator '%s'"%self.operator) 161 return logOp.join( 162 '(%s)'%self._flatLogic(template, op) for op in self.operands)

163

164 165 -def _make_pgsposition(coords):

166 if len(coords)!=2: 167 raise common.STCSParseError("STC-S points want two coordinates.") 168 return pgsphere.SPoint(*coords)

169

170 171 -def _make_pgscircle(coords):

172 if len(coords)!=3: 173 raise common.STCSParseError("STC-S circles want three numbers.") 174 return pgsphere.SCircle(pgsphere.SPoint(*coords[:2]), coords[2])

175

176 177 -def _make_pgsbox(coords):

178 if len(coords)!=4: 179 raise common.STCSParseError("STC-S boxes want four numbers.") 180 x,y,w,h = coords 181 return pgsphere.SPoly(( 182 pgsphere.SPoint(x-w/2, y-h/2), 183 pgsphere.SPoint(x-w/2, y+h/2), 184 pgsphere.SPoint(x+w/2, y+h/2), 185 pgsphere.SPoint(x+w/2, y-h/2)))

186

187 188 -def _make_pgspolygon(coords):

189 if len(coords)<6 or len(coords)%2: 190 raise common.STCSParseError( 191 "STC-S polygons want at least three number pairs") 192 return pgsphere.SPoly( 193 [pgsphere.SPoint(*p) for p in utils.iterConsecutivePairs(coords)])

194

195 196 -def _makePgSphereInstance(match):

197 """returns a utils.pgsphere instance from a match of simpleStatement in 198 the simple STCS parser below. 199 """ 200 if match["flavor"] and match["flavor"].strip().upper()!="SPHERICAL2": 201 raise common.STCSParseError("Only SPHERICAL2 STC-S supported here") 202 refFrame = 'UnknownFrame' 203 if match["frame"]: 204 refFrame = match["frame"].strip() 205 # refFrame gets thrown away here; to use it, we'd have to generate 206 # ADQL nodes, and that would be clumsy for uploads. See rant above. 207 handler = globals()["_make_pgs%s"%match["shape"].lower()] 208 res = handler( 209 tuple(float(s)*utils.DEG for s in match["coords"].strip().split() if s)) 210 res.cooSys = refFrame 211 return res

212

213 -def _makeRE(strings):

214 return "|".join(sorted(strings, key=lambda s: -len(s)))

215

216 @utils.memoized 217 -def getSimpleSTCSParser():

218 from pyparsing import ( 219 Regex, CaselessKeyword, OneOrMore, Forward, Suppress, 220 Optional, ParseException, ParseSyntaxException) 221 222 with utils.pyparsingWhitechars(" \t\n\r"): 223 frameRE = _makeRE(TAP_SYSTEMS) 224 refposRE = _makeRE(TAP_REFPOS) 225 flavorRE = _makeRE(TAP_FLAVORS) 226 systemRE = (r"(?i)\s*" 227 r"(?P<frame>%s)?\s*" 228 r"(?P<refpos>%s)?\s*" 229 r"(?P<flavor>%s)?\s*")%( 230 frameRE, refposRE, flavorRE) 231 coordsRE = r"(?P<coords>(%s\s*)+)"%utils.floatRE 232 233 simpleStatement = Regex("(?i)\s*" 234 "(?P<shape>position|circle|box|polygon)" 235 +systemRE 236 +coordsRE) 237 simpleStatement.setName("STC-S geometry") 238 simpleStatement.addParseAction(lambda s,p,t: _makePgSphereInstance(t)) 239 system = Regex(systemRE) 240 system.setName("STC-S system spec") 241 region = Forward() 242 notExpr = CaselessKeyword("NOT") + Suppress('(') + region + Suppress(')') 243 notExpr.addParseAction(lambda s,p,t: GeomExpr("UNKNOWN", "NOT", (t[1],))) 244 opExpr = ( 245 (CaselessKeyword("UNION") | CaselessKeyword("INTERSECTION"))("op") 246 + Optional(Regex(frameRE))("frame") 247 + Optional(Regex(refposRE)) + Optional(Regex(flavorRE)) 248 + Suppress("(") 249 + region + OneOrMore(region) 250 + Suppress(")")) 251 opExpr.addParseAction( 252 lambda s,p,t: GeomExpr(str(t["frame"]), t[0].upper(), t[2:])) 253 region << (simpleStatement | opExpr | notExpr) 254 255 def parse(s): 256 if isinstance(s, pgsphere.PgSAdapter): 257 return s 258 if s is None or not s.strip(): # special service: Null values 259 return None 260 261 try: 262 res = utils.pyparseString(region, s, parseAll=True)[0] 263 if not res.cooSys or res.cooSys.lower()=='unknownframe': # Sigh. 264 res.cooSys = "UNKNOWN" 265 return res 266 except (ParseException, ParseSyntaxException) as msg: 267 raise common.STCSParseError("Invalid STCS (%s)"%str(msg))

268 269 return parse 270 271 272 parseSimpleSTCS = getSimpleSTCSParser()

273 274 275 -def simpleSTCSToPolygon(stcsExpr):

276 """returns an spoly for an STC-S region specification. 277 278 This is used to let people upload VOTables with REGION items. 279 """ 280 if stcsExpr is None: 281 return None 282 return parseSimpleSTCS(stcsExpr).asPoly()

283 284 285 if __name__=="__main__": 286 # compute the Euler angles given above. pgsphere has its rotation 287 # matrices opposite to ours (ccw rather than cw), hence the negation 288 from kapteyn import celestial 289 from gavo.stc import sphermath 290 print("FK4:", tuple(-a for a in 291 sphermath.getEulerAnglesFromMatrix( 292 celestial.skymatrix("icrs", "fk4 B1950.0")[0]))) 293 print("GAL:", tuple(-a for a in 294 sphermath.getEulerAnglesFromMatrix( 295 celestial.skymatrix("icrs", "galactic")[0]))) 296

Source Code for Module gavo.stc.tapstc