gavo.stc.bboxes

1 """ 2 Computing bboxes for STC geometries. 3 4 A bbox coming out of this module is a 4-tuple of (ra0, de0, ra1, de1) in 5 ICRS degrees. 6 7 (You're right, this should be part of the dm classes; but it's enough 8 messy custom code that I found it nicer to break it out). 9 """ 10 11 #c Copyright 2008-2019, the GAVO project 12 #c 13 #c This program is free software, covered by the GNU GPL. See the 14 #c COPYING file in the source distribution. 15 16 17 import math 18 from math import sin, cos, tan, atan, acos 19 20 from gavo import utils 21 from gavo.stc import common 22 from gavo.stc import conform 23 from gavo.stc import stcsast 24 25 from gavo.utils import DEG

26 27 ############### Some spherical geometry calculations 28 # (most due to Chamberlain and Duquette, "Some Algorithms for Polygons on 29 # a sphere", http://hdl.handle.net/2014/40409) 30 31 32 -def getHeading(long1, lat1, long2, lat2):

33 """returns the initial heading when going from one spherical position to 34 another along a great circle. 35 36 Everything is in rad, angle is counted north over east. 37 """ 38 if abs(long1-long2)<=1e-10: 39 # along a meridian 40 if lat1<lat2: 41 return 0 42 else: 43 return math.pi 44 45 if lat1<=-math.pi/2+1e-10: 46 # starting on south pole 47 return 0 48 elif lat1>=math.pi/2-1e-10: 49 # starting on north pole 50 return math.pi 51 52 zeta = utils.spherDist( 53 utils.spherToCart(long1, lat1), 54 utils.spherToCart(long2, lat2)) 55 cosPhi = (sin(lat2)-sin(lat1)*cos(zeta))/(sin(zeta)*cos(lat1)) 56 if abs(cosPhi-1)<1e-10: 57 return 0 58 elif abs(cosPhi+1)<1e-10: 59 return math.pi 60 elif sin(long2-long1)>0: 61 return acos(cosPhi) 62 else: 63 return common.TWO_PI-acos(cosPhi)

64

65 66 -class GCSegment(object):

67 """A great circle segment on a sphere. It is assumed that it's no larger 68 than pi. 69 70 Construction is long,lat, long, lat in rad (or use the fromDegrees 71 class method. 72 73 Very small (<1e-8 rad, say) segments are not supported. You should 74 probably work in the tangential plane on that kind of scale. 75 """

76 - def __init__(self, long1, lat1, long2, lat2):

77 self.long1, self.lat1 = long1, lat1 78 self.long2, self.lat2 = long2, lat2 79 self._normalize() 80 if self.long1>self.long2: 81 self.long1, self.long2 = self.long2, self.long1 82 self.lat1, self.lat2 = self.lat2, self.lat1 83 84 vertical = abs(self.long1-self.long2)<1e-10 85 if vertical and abs(self.lat1-self.lat2)<1e-10: 86 raise ValueError("Null segment: start and end are identical") 87 88 self.overPole = vertical or cos(lat1)*cos(lat2)<1e-10 89 self.overStich = 2*math.pi<=self.long1+self.long2<3*math.pi

90 91 @classmethod

92 - def fromDegrees(cls, long1, lat1, long2, lat2):

93 return cls(long1*DEG, lat1*DEG, long2*DEG, lat2*DEG)

94

95 - def __str__(self):

96 return "<Great Circle through (%f %f), (%f %f)"%( 97 self.long1/DEG, self.lat1/DEG, self.long2/DEG, self.lat2/DEG)

98

99 - def _normalize(self):

100 """makes sure all longitudes are in [0, 2 pi[ and all latitudes between 101 [-pi, pi]. 102 """ 103 self.long1 = common.clampLong(self.long1) 104 self.long2 = common.clampLong(self.long2) 105 self.lat1 = common.clampLat(self.lat1) 106 self.lat2 = common.clampLat(self.lat2)

107

108 - def _computeBBFor(self, long1, lat1, long2, lat2):

109 """returns a 4-tuple bounding box for a great circle segment. 110 111 This is restricted to GCs not crossing the stitch line; furthermore, it 112 still depends on self. 113 """ 114 if self.overPole: 115 latMin, latMax = min(lat1, lat2), max(lat1, lat2) 116 if cos(lat1)*cos(lat2)<1e-10: 117 # the pole is part of the segment 118 long1, long2 = 0, common.TWO_PI 119 120 else: 121 # compute headings at both ends 122 theta1 = getHeading(long1, lat1, long2, lat2) 123 theta2 = getHeading(long2, lat2, long1, lat1) 124 125 # if the headings are two quadrants apart, the limits are at 126 # the corners 127 quad1, quad2 = int(2*theta1/math.pi), int(2*theta2/math.pi) 128 if abs(quad1-quad2)==2: 129 latExt = lat1 130 else: 131 # the limit is somewhere in between. 132 latExt = acos(abs(sin(theta1))*cos(lat1)) 133 # acos has two branches; inspecting lat1 for the sign is 134 # enough since latExt won't be extreme when sgn(lat1)!=sgn(lat2). 135 if lat1<0: 136 latExt = -latExt 137 138 latMin = min(latExt, lat1, lat2) 139 latMax = max(latExt, lat1, lat2) 140 return (long1, latMin, long2, latMax)

141

142 - def latForLong(self, long):

143 """returns the latitude the great circle is on for a given longitude. 144 145 Input and output is in rad. 146 147 If the pole is part of the great circle, a constant (but fairly 148 meaningless) value is returned. 149 """ 150 if self.overPole: 151 return self.long1 152 return atan( 153 tan(self.lat1)*(sin(long-self.long2)/sin(self.long1-self.long2)) 154 -tan(self.lat2)*(sin(long-self.long1)/sin(self.long1-self.long2)))

155

156 - def getBBs(self):

157 """returns a sequence of bounding boxes for this great circle segment. 158 159 Actually, the sequence contains one box for a segment that does not 160 cross the stitching line, two boxes for those that do. 161 """ 162 if not self.overPole and common.TWO_PI<=self.long1+self.long2<3*math.pi: 163 # segment crosses stitch line, split it up 164 latAtStich = self.latForLong(0) 165 return [ 166 self._computeBBFor(0, latAtStich, self.long1, self.lat1), 167 self._computeBBFor(self.long2, self.lat2, 2*math.pi, latAtStich)] 168 else: 169 return [ 170 self._computeBBFor(self.long1, self.lat1, self.long2, self.lat2)]

171

172 173 ############################################################# 174 # Actual bbox hacks. I'm pretty sure there's loads of bad 175 # corner cases in here. In particular in the vicinity of the 176 # pole, *much* more thought and testing is required (I'm hoping 177 # astropy will come to the rescue at some point). 178 179 @utils.memoized 180 -def getStandardFrame():

181 return stcsast.parseSTCS("Position ICRS unit deg")

182

183 184 -def _makeSphericalBbox(minRA, minDec, maxRA, maxDec):

185 """yields one or two bboxes from for spherical coordinates. 186 187 This handles crossing the stich line as well as shooting over the pole. 188 189 Everything here is in degrees. 190 191 This function assumes that -360<minRA<maxRA<720 and that at least one 192 of minRA, maxRA is between 0 and 360. 193 """ 194 if minDec<-90: 195 # fold over the pole 196 maxDec = max(maxDec, -minDec-180) 197 minDec = -90 198 minRA, maxRA = 0, 360 199 200 if maxDec>90: 201 # fold over the pole 202 minDec = min(minDec, 180-maxDec) 203 maxDec = 90 204 minRA, maxRA = 0, 360 205 206 if minRA<0: 207 yield (0, minDec, maxRA, maxDec) 208 yield (360+minRA, minDec, 360, maxDec) 209 elif maxRA>360: 210 yield (0, minDec, maxRA-360, maxDec) 211 yield (minRA, minDec, 360, maxDec) 212 else: 213 yield (minRA, minDec, maxRA, maxDec)

214

215 216 -def _intersectBboxes(bbox1, bbox2):

217 """returns the intersection of the two bboxes. 218 """ 219 return ( 220 max(bbox1[0], bbox2[0]), 221 max(bbox1[1], bbox2[1]), 222 min(bbox1[2], bbox2[2]), 223 min(bbox1[3], bbox2[3]))

224

225 226 -def _isEmpty(bbox):

227 """returns true if bbox is empty. 228 """ 229 return bbox[0]>=bbox[2] or bbox[1]>=bbox[3]

230

231 232 -def _computeIntersectionForSequence(seq1, seq2):

233 """helps _computeIntersection; see comment there. 234 """ 235 if seq1 is None: 236 return seq2 237 if seq2 is None: 238 return seq1 239 240 commonBoxes = [] 241 for box1 in seq1: 242 for box2 in seq2: 243 commonBox = _intersectBboxes(box1, box2) 244 if not _isEmpty(commonBox): 245 commonBoxes.append(commonBox) 246 return commonBoxes

247

248 249 -def _computeCircleBbox(circle):

250 """helps _getBboxesFor. 251 """ 252 return _makeSphericalBbox( 253 circle.center[0]-circle.radius, 254 circle.center[1]-circle.radius, 255 circle.center[0]+circle.radius, 256 circle.center[1]+circle.radius)

257

258 259 -def _computeEllipseBbox(ellipse):

260 """helps _getBboxesFor. 261 262 TODO: Make tighter by actually using minor axis. 263 """ 264 return _makeSphericalBbox( 265 ellipse.center[0]-ellipse.smajAxis, 266 ellipse.center[1]-ellipse.smajAxis, 267 ellipse.center[0]+ellipse.smajAxis, 268 ellipse.center[1]+ellipse.smajAxis)

269

270 271 -def _computeBoxBbox(box):

272 """helps _getBboxesFor. 273 274 Note that the box boundaries are great circles. 275 """ 276 maxDec = max(-90, min(90, box.center[1]+box.boxsize[1])) 277 minDec = max(-90, min(90, box.center[1]-box.boxsize[1])) 278 279 if maxDec>90-1e-10: 280 yield (0, minDec, 360, 90) 281 return 282 if minDec<-90+1e-10: 283 yield (0, -90, 360, maxDec) 284 return 285 286 upperBBs = GCSegment.fromDegrees( 287 box.center[0]-box.boxsize[0], maxDec, 288 box.center[0]+box.boxsize[0], maxDec 289 ).getBBs() 290 lowerBBs = GCSegment.fromDegrees( 291 box.center[0]-box.boxsize[0], minDec, 292 box.center[0]+box.boxsize[0], minDec 293 ).getBBs() 294 295 for bbox in _makeSphericalBbox( 296 *joinBboxes(fromRad(upperBBs[0]), fromRad(lowerBBs[0]))): 297 yield bbox 298 if len(upperBBs)==2: # stitching required 299 assert len(lowerBBs)==2 300 for bbox in _makeSphericalBbox( 301 *joinBboxes(fromRad(upperBBs[1]), fromRad(lowerBBs[1]))): 302 yield bbox

303

304 305 -def _computePolygonBbox(poly):

306 """helps _getBboxesFor. 307 """ 308 leftBoxes, rightBoxes = [], [] 309 310 def addBoxes(b): 311 rightBoxes.append(fromRad(b[0])) 312 if len(b)>1: 313 leftBoxes.append(fromRad(b[1]))

314 315 firstPoint = lastPoint = poly.vertices[0] 316 for point in poly.vertices[1:]: 317 curSegment = GCSegment.fromDegrees( 318 lastPoint[0], lastPoint[1], point[0], point[1]) 319 addBoxes(list(curSegment.getBBs())) 320 lastPoint = point 321 322 try: 323 curSegment = GCSegment.fromDegrees( 324 lastPoint[0], lastPoint[1], firstPoint[0], firstPoint[1]) 325 except ValueError: # probably a null segment since the authors 326 # closed the polygon. Ignore this. 327 pass 328 else: 329 addBoxes(list(curSegment.getBBs())) 330 331 yield joinBboxes(*rightBoxes) 332 if leftBoxes: 333 yield joinBboxes(*leftBoxes) 334

335 336 -def _computeAllSkyBbox(geo):

337 """helps _getBboxesFor. 338 """ 339 yield (0, -90, 360, 90)

340

341 342 -def _computeSpaceIntervalBbox(geo):

343 """helps _getBboxesFor. 344 345 PositionInterval could have all kinds of weird coordinate systems. 346 We only check there's exactly two dimensions and otherwise hope for the 347 best. 348 """ 349 if geo.frame.nDim!=2: 350 raise ValueError("Only PositionsIntervals with nDim=2 are supported" 351 " for bboxes.") 352 353 lowerLimit = geo.lowerLimit 354 if lowerLimit is None: 355 lowerLimit = (0, -90) 356 357 upperLimit = geo.upperLimit 358 if upperLimit is None: 359 upperLimit = (360, 90) 360 return _makeSphericalBbox(lowerLimit[0], lowerLimit[1], 361 upperLimit[0], upperLimit[1])

362

363 364 -def _computeUnionBbox(geo):

365 """helps _getBboxesFor. 366 367 Union is just returning all bboxes from our child geometries. 368 """ 369 for child in geo.children: 370 for bbox in _getBboxesFor(child): 371 yield bbox

372

373 374 -def _computeIntersectionBbox(geo):

375 """helps _getBboxesFor. 376 """ 377 # This is surprisingly involved, since our children may yield an arbitrary 378 # number of part boxes; each of those could contribute to an intersection. 379 # To figure out the intersection, we need to intersect each part bbox with 380 # each other part bbox and keep whatever is non-empty. 381 commonBboxes = None 382 for child in geo.children: 383 commonBboxes = _computeIntersectionForSequence(commonBboxes, 384 list(_getBboxesFor(child))) 385 for bbox in commonBboxes: 386 yield bbox

387

388 -def _computeDifferenceBbox(geo):

389 """helps _getBboxesFor. 390 391 (we ignore the cut-out). 392 """ 393 return _getBboxesFor(geo.children[0])

394

395 -def _computeNotBbox(geo):

396 """helps _getBboxesFor. 397 398 (we ignore the cut-out and always return the entire sky) 399 """ 400 # of course, we could provide a saner implementation of _computeIntersection 401 # which then would make up to eight bboxes out of such a thing. 402 # Never mind, people shouldn't do this anyway. 403 return _computeAllSkyBbox(geo)

404 405 406 407 _BBOX_COMPUTERS = { 408 "Circle": _computeCircleBbox, 409 "Ellipse": _computeEllipseBbox, 410 "Box": _computeBoxBbox, 411 "Polygon": _computePolygonBbox, 412 "AllSky": _computeAllSkyBbox, 413 "SpaceInterval": _computeSpaceIntervalBbox, 414 "Union": _computeUnionBbox, 415 "Intersection": _computeIntersectionBbox, 416 "Difference": _computeDifferenceBbox, 417 "Not": _computeNotBbox, 418 }

419 420 421 -def _getBboxesFor(geo):

422 """yields one or two bboxes for a conformed geometry. 423 424 Two bboxes are returned when the geometry overlaps the stitching point. 425 426 A geometry is conformed if it's been conformed to what's coming back 427 from getStandardFrame() above. 428 """ 429 geoName = geo.__class__.__name__ 430 if geoName not in _BBOX_COMPUTERS: 431 raise common.STCInternalError("Do now know how to compute the bbox of" 432 " %s."%geoName) 433 434 for bbox in _BBOX_COMPUTERS[geoName](geo): 435 yield bbox

436

437 438 -def fromRad(bbox):

439 return tuple(a/utils.DEG for a in bbox)

440

441 442 -def joinBboxes(*bboxes):

443 """returns a bounding box encompassing all the the bboxes passed in. 444 445 No input bbox must cross the stitching line; they must be normalized, 446 i.e., lower left and upper right corners. 447 """ 448 if not bboxes: 449 raise common.InternalError("bbox join without bbox") 450 minRA, maxRA = utils.Supremum, utils.Infimum 451 minDec, maxDec = utils.Supremum, utils.Infimum 452 453 for ra0, dec0, ra1, dec1 in bboxes: 454 minRA = min(minRA, ra0) 455 minDec = min(minDec, dec0) 456 maxRA = max(maxRA, ra1) 457 maxDec = max(maxDec, dec1) 458 459 return (minRA, minDec, maxRA, maxDec)

460

461 462 -def getBboxes(ast):

463 """iterates over the bboxes of the areas within ast. 464 465 bboxes are (ra0, de0, ra1, de1) in ICRS degrees. 466 """ 467 astc = conform.conform(ast, getStandardFrame()) 468 for area in astc.areas: 469 for bbox in _getBboxesFor(area): 470 yield bbox

471

Source Code for Module gavo.stc.bboxes