gavo.protocols.siap

1 """ 2 Support code for the Simple Image Access Protocol. 3 """ 4 5 #c Copyright 2008-2019, the GAVO project 6 #c 7 #c This program is free software, covered by the GNU GPL. See the 8 #c COPYING file in the source distribution. 9 10 11 import math 12 import re 13 import urllib 14 15 from gavo import base 16 from gavo import svcs 17 from gavo import utils 18 from gavo.base import coords 19 from gavo.protocols import products 20 from gavo.utils import fitstools 21 from gavo.utils import pgsphere 22 23 MS = base.makeStruct 24 25 ####################### SIAPv2 magic 26

27 -def parseSIAP2Geometry(aString, fieldName="POS"):

28 """parses a SIAPv2 geometry spec to a pgsphere object. 29 30 Parse errors raise validation errors for fieldName. 31 """ 32 mat = re.match("(CIRCLE|RANGE|POLYGON) (.*)", aString) 33 if not mat: 34 raise base.ValidationError("Invalid SIAPv2 geometry: '%s'" 35 " (expected a SIAPv2 shape name)"%utils.makeEllipsis(aString, 20), 36 fieldName) 37 38 geoName = mat.group(1) 39 try: 40 args = [float(s) for s in mat.group(2).split()] 41 except ValueError: 42 raise base.ValidationError("Invalid SIAPv2 coordinates: '%s'" 43 " (bad floating point literal '%s')"%( 44 utils.makeEllipsis(mat.group(2), 20), s), 45 fieldName) 46 47 if geoName=="CIRCLE": 48 if len(args)!=3: 49 raise base.ValidationError("Invalid SIAPv2 CIRCLE: '%s'" 50 " (need exactly three numbers)"%( 51 utils.makeEllipsis(aString, 20)), 52 fieldName) 53 return pgsphere.SCircle(pgsphere.SPoint.fromDegrees(args[0], args[1]), 54 args[2]*utils.DEG) 55 56 elif geoName=="RANGE": 57 # SBox isn't really RANGE, but RANGE shouldn't have been 58 # part of the standard and people that use it deserve 59 # to get bad results. 60 if len(args)!=4: 61 raise base.ValidationError("Invalid SIAPv2 RANGE: '%s'" 62 " (need exactly four numbers)"%( 63 utils.makeEllipsis(aString, 20)), 64 fieldName) 65 if args[0]>args[1] or args[2]>args[3]: 66 raise base.ValidationError("Invalid SIAPv2 RANGE: '%s'" 67 " (lower limits must be smaller than upper limits)"%( 68 utils.makeEllipsis(aString, 20)), 69 fieldName) 70 return pgsphere.SBox( 71 pgsphere.SPoint.fromDegrees(args[0], args[2]), 72 pgsphere.SPoint.fromDegrees(args[1], args[3])) 73 74 elif geoName=="POLYGON": 75 if len(args)<6 or len(args)%2: 76 raise base.ValidationError("Invalid SIAPv2 POLYGON: '%s'" 77 " (need more than three coordinate *pairs*)"%( 78 utils.makeEllipsis(mat.group(2), 20)), 79 fieldName) 80 return pgsphere.SPoly([ 81 pgsphere.SPoint.fromDegrees(*pair) 82 for pair in utils.iterConsecutivePairs(args)]) 83 84 else: 85 assert False

86 87 88 ####################### pgsSIAP mixin 89 90 # expressions as used in getPGSQuery 91 _PGS_OPERATORS = { 92 "COVERS": "coverage ~ %%(%s)s", 93 "ENCLOSED": "%%(%s)s ~ coverage", 94 "CENTER": None, # special handling below 95 "OVERLAPS": "%%(%s)s && coverage", 96 } 97

98 -def getPGSQuery(intersect, ra, dec, sizes, prefix, sqlPars):

99 """returns SQL for a SIAP query on pgsSIAP tables. 100 """ 101 if intersect=='CENTER': 102 return "%%(%s)s @ coverage"%(base.getSQLKey( 103 prefix+"center", pgsphere.SPoint.fromDegrees(ra, dec), sqlPars)) 104 105 expr = _PGS_OPERATORS[intersect] 106 try: 107 targetBox = pgsphere.SBox.fromSIAPPars(ra, dec, sizes[0], sizes[1]) 108 return expr%base.getSQLKey(prefix+"area", targetBox, sqlPars) 109 except pgsphere.TwoSBoxes as ex: 110 # Fold-over at pole, return a disjunction 111 return "( %s OR %s )"%( 112 expr%base.getSQLKey(prefix+"area1", ex.box1, sqlPars), 113 expr%base.getSQLKey(prefix+"area2", ex.box2, sqlPars))

114 115 116 ####################### SIAP service helpers, cores, etc. 117

118 -def dissectPositions(posStr):

119 """tries to infer RA and DEC from posStr. 120 121 In contrast to base.parseCooPair, we are quite strict here and just 122 try to cope with some bad clients that leave out the comma. 123 """ 124 try: 125 ra, dec = map(float, posStr.split(",")) 126 except ValueError: # maybe a sign as separator? 127 if '+' in posStr: 128 ra, dec = map(float, posStr.split("+")) 129 elif '-' in posStr: 130 ra, dec = map(float, posStr.split("-")) 131 else: 132 raise ValueError("No pos") 133 return ra, dec

134 135

136 -def _getQueryMaker(queriedTable):

137 """returns a query making function for SIAP appropriate for queriedTable. 138 139 This used to have a function when we had different backends for SIAP. 140 Curently, we no longer have that, so this always returns getPGSQuery. 141 """ 142 return getPGSQuery

143 144

145 -def getQuery(queriedTable, parameters, sqlPars, prefix="sia"):

146 """returns an SQL fragment for a SIAP query for bboxes. 147 148 The SQL is returned as a WHERE-fragment in a string. The parameters 149 are added in the sqlPars dictionary. 150 151 parameters is a dictionary that maps the SIAP keywords to the 152 values in the query. Parameters not defined by SIAP are ignored. 153 """ 154 posStr = urllib.unquote(parameters["POS"]) 155 try: 156 ra, dec = dissectPositions(posStr) 157 except (ValueError, TypeError): 158 raise base.ui.logOldExc(base.ValidationError( 159 "%s is not a RA,DEC pair."%posStr, "POS", posStr)) 160 try: 161 sizes = map(float, parameters["SIZE"].split(",")) 162 except ValueError: 163 raise base.ui.logOldExc(base.ValidationError("Size specification" 164 " has to be <degs> or <degs>,<degs>", "SIZE", parameters["SIZE"])) 165 if len(sizes)==1: 166 sizes = sizes*2 167 intersect = parameters.get("INTERSECT", "OVERLAPS") 168 query = _getQueryMaker(queriedTable)( 169 intersect, ra, dec, sizes, prefix, sqlPars) 170 # the following are for the benefit of cutout queries. 171 sqlPars["_ra"], sqlPars["_dec"] = ra, dec 172 sqlPars["_sra"], sqlPars["_sdec"] = sizes 173 return query

174 175

176 -class SIAPCutoutCore(svcs.DBCore):

177 """A core doing SIAP plus cutouts. 178 179 It has, by default, an additional column specifying the desired size of 180 the image to be retrieved. Based on this, the cutout core will tweak 181 its output table such that references to cutout images will be retrieved. 182 183 The actual process of cutting out is performed by the product core and 184 renderer. 185 """ 186 name_ = "siapCutoutCore" 187 188 # This should become a property or something once we 189 # compress the stuff or have images with bytes per pixel != 2 190 copiedCols = ["centerAlpha", "centerDelta", "imageTitle", "instId", 191 "dateObs", "nAxes", "pixelSize", "pixelScale", "mime", 192 "refFrame", "wcs_equinox", "wcs_projection", "wcs_refPixel", 193 "wcs_refValues", "wcs_cdmatrix", "bandpassId", "bandpassUnit", 194 "bandpassHi", "bandpassLo", "pixflags"] 195

196 - def getQueryCols(self, service, queryMeta):

197 cols = svcs.DBCore.getQueryCols(self, service, queryMeta) 198 for name in self.copiedCols: 199 cols.append(svcs.OutputField.fromColumn( 200 self.queriedTable.getColumnByName(name))) 201 d = self.queriedTable.getColumnByName("accsize").copy(self) 202 d.tablehead = "Est. file size" 203 cols.append(svcs.OutputField.fromColumn(d)) 204 return cols

205

206 - def _fixRecord(self, record, centerAlpha, centerDelta, sizeAlpha, sizeDelta):

207 """inserts estimates for WCS values into a cutout record. 208 """ 209 origWCS = coords.getWCS({ 210 "CUNIT1": "deg", "CUNIT2": "deg", "CTYPE1": "RA---TAN", 211 "CTYPE2": "DEC--TAN", 212 "CRVAL1": record["wcs_refValues"][0], 213 "CRVAL2": record["wcs_refValues"][1], 214 "CRPIX1": record["wcs_refPixel"][0], 215 "CRPIX2": record["wcs_refPixel"][1], 216 "CD1_1": record["wcs_cdmatrix"][0], 217 "CD1_2": record["wcs_cdmatrix"][1], 218 "CD2_1": record["wcs_cdmatrix"][2], 219 "CD2_2": record["wcs_cdmatrix"][3], 220 "LONPOLE": 180., 221 "NAXIS": record["nAxes"], 222 "NAXIS1": record["pixelSize"][0], 223 "NAXIS2": record["pixelSize"][1], 224 }) 225 limits = coords.getPixelLimits( 226 ((centerAlpha-sizeAlpha/2, centerDelta-sizeDelta/2), 227 (centerAlpha+sizeAlpha/2, centerDelta+sizeDelta/2)), 228 origWCS) 229 _, newHeader = fitstools.cutoutHeader(origWCS._dachs_header, *limits) 230 231 if newHeader["NAXIS1"]<3 or newHeader["NAXIS2"]<3: 232 # these typically result from crazy calibrations or very marginal 233 # matches that nobody wants to see 234 raise utils.SkipThis("Too marginal a match") 235 236 record["wcs_refPixel"] = [newHeader["CRPIX1"], newHeader["CRPIX2"]] 237 record["wcs_refValues"] = [newHeader["CRVAL1"], newHeader["CRVAL2"]] 238 record["pixelSize"] = [newHeader["NAXIS1"], newHeader["NAXIS2"]] 239 240 record["accref"] = products.RAccref(record["accref"], { 241 "ra": centerAlpha, "dec": centerDelta, 242 "sra": sizeAlpha, "sdec": sizeDelta}) 243 record["centerAlpha"] = centerAlpha 244 record["centerDelta"] = centerDelta 245 246 newWCS = coords.getWCS(newHeader) 247 corners = newWCS.calcFootprint() 248 record["coverage"] = pgsphere.SPoly.fromDALI([ 249 corners[0][0], corners[0][1], 250 corners[1][0], corners[1][1], 251 corners[2][0], corners[2][1], 252 corners[3][0], corners[3][1]]) 253 # Aw, fix this already. 254 bytesPerPixel = 2 255 record["accsize"] = min(record["accsize"], 256 int(bytesPerPixel*newHeader["NAXIS1"]*newHeader["NAXIS2"]))

257

258 - def run(self, service, inputData, queryMeta):

259 res = svcs.DBCore.run(self, service, inputData, queryMeta) 260 sqlPars = queryMeta["sqlQueryPars"] 261 try: 262 sra = sdec = float(queryMeta.ctxArgs["cutoutSize"]) 263 except (KeyError, ValueError): 264 try: 265 sra, sdec = sqlPars["_sra"], sqlPars["_sdec"] 266 except KeyError: 267 sra, sdec = 0.5, 0.5 268 269 if "_dec" in sqlPars: 270 cosD = math.cos(sqlPars["_dec"]/180*math.pi) 271 if abs(cosD)>1e-5: 272 sra = sra/cosD 273 else: 274 sra = 360 275 276 newRows = [] 277 for record in res: 278 try: 279 self._fixRecord(record, 280 sqlPars.get("_ra", record["centerAlpha"]), 281 sqlPars.get("_dec", record["centerDelta"]), sra, sdec) 282 newRows.append(record) 283 except base.SkipThis: 284 pass 285 except ValueError as msg: 286 # Old pywcs signified troublesome WCS with a ValueError-derived 287 # thing. Let's see how this works out for astropy. 288 base.ui.notifyWarning("Botched WCS (%s) in the record %s"%( 289 msg, record)) 290 291 res.rows = newRows 292 293 return res

294 295

296 -def _test():

297 import doctest, siap 298 doctest.testmod(siap)

299 300 301 if __name__=="__main__": 302 _test() 303

Source Code for Module gavo.protocols.siap