Coverage for python/lsst/skymap/tractInfo.py: 30%

2# LSST Data Management System

5# This product includes software developed by the

6# LSST Project (http://www.lsst.org/).

8# This program is free software: you can redistribute it and/or modify

9# it under the terms of the GNU General Public License as published by

10# the Free Software Foundation, either version 3 of the License, or

11# (at your option) any later version.

12#

13# This program is distributed in the hope that it will be useful,

14# but WITHOUT ANY WARRANTY; without even the implied warranty of

15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

16# GNU General Public License for more details.

17#

18# You should have received a copy of the LSST License Statement and

19# the GNU General Public License along with this program. If not,

20# see <http://www.lsstcorp.org/LegalNotices/>.

21#

23__all__ = ["TractInfo"]

25import numpy as np

26from deprecated.sphinx import deprecated

28import lsst.pex.exceptions

29import lsst.geom as geom

30from lsst.sphgeom import ConvexPolygon, Box

32from .detail import makeSkyPolygonFromBBox, Index2D

35class TractInfo:

36 """Information about a tract in a SkyMap sky pixelization

38 Parameters

39 ----------

40 id : `int`

41 tract ID

42 tractBuilder : Subclass of `lsst.skymap.BaseTractBuilder`

43 Object used to compute patch geometry.

44 ctrCoord : `lsst.geom.SpherePoint`

45 ICRS sky coordinate of center of inner region of tract; also used as

46 the CRVAL for the WCS.

47 vertexCoordList : `list` of `lsst.geom.SpherePoint`

48 Vertices that define the boundaries of the inner region.

49 tractOverlap : `lsst.geom.Angle`

50 Minimum overlap between adjacent sky tracts; this defines the minimum

51 distance the tract extends beyond the inner region in all directions.

52 wcs : `lsst.afw.image.SkyWcs`

53 WCS for tract. The reference pixel will be shifted as required so that

54 the lower left-hand pixel (index 0,0) has pixel position 0.0, 0.0.

55 innerBoxCorners : `list` [`lsst.sphgeom.LonLat`], optional

56 If set then the ``inner_sky_region`` will be a `lsst.sphgeom.Box` with

57 these corners as opposed to a `lsst.sphgeom.ConvexPolygon` built from

58 the ``vertex_list``.

60 Notes

61 -----

62 The tract is subdivided into rectangular patches. Each patch has the

63 following properties:

65 - An inner region defined by an inner bounding box. The inner regions of

66 the patches exactly tile the tract, and all inner regions have the same

67 dimensions. The tract is made larger as required to make this work.

69 - An outer region defined by an outer bounding box. The outer region

70 extends beyond the inner region by patchBorder pixels in all directions,

71 except there is no border at the edges of the tract.

72 Thus patches overlap each other but never extend off the tract.

73 If you do not want any overlap between adjacent patches then set

74 patchBorder to 0.

76 - An index that consists of a pair of integers:

78 * 0 <= x index < numPatches[0]

80 * 0 <= y index < numPatches[1]

82 Patch 0,0 is at the minimum corner of the tract bounding box.

84 - It is not enforced that ctrCoord is the center of vertexCoordList, but

85 SkyMap relies on it.

86 """

87 def __init__(self, id, tractBuilder, ctrCoord, vertexCoordList, tractOverlap, wcs, innerBoxCorners=None):

88 self._id = id

89 self._ctrCoord = ctrCoord

90 self._vertexCoordList = tuple(vertexCoordList)

91 self._tractOverlap = tractOverlap

92 self._tractBuilder = tractBuilder

93 self._innerBoxCorners = innerBoxCorners

95 minBBox = self._minimumBoundingBox(wcs)

96 initialBBox, self._numPatches = self._tractBuilder.setupPatches(minBBox, wcs)

97 self._bbox, self._wcs = self._finalOrientation(initialBBox, wcs)

99 def _minimumBoundingBox(self, wcs):

100 """Calculate the minimum bounding box for the tract, given the WCS.

101

102 The bounding box is created in the frame of the supplied WCS,

103 so that it's OK if the coordinates are negative.

104

105 We compute the bounding box that holds all the vertices and the

106 desired overlap.

107 """

108 minBBoxD = geom.Box2D()

109 halfOverlap = self._tractOverlap / 2.0

110 for vertexCoord in self._vertexCoordList:

111 if self._tractOverlap == 0:

112 minBBoxD.include(wcs.skyToPixel(vertexCoord))

113 else:

114 numAngles = 24

115 angleIncr = geom.Angle(360.0, geom.degrees) / float(numAngles)

116 for i in range(numAngles):

117 offAngle = angleIncr * i

118 offCoord = vertexCoord.offset(offAngle, halfOverlap)

119 pixPos = wcs.skyToPixel(offCoord)

120 minBBoxD.include(pixPos)

121 return minBBoxD

122

123 def _finalOrientation(self, bbox, wcs):

124 """Determine the final orientation

125

126 We offset everything so the lower-left corner is at 0,0

127 and compute the final Wcs.

128

129 Parameters

130 ----------

131 bbox : `lsst.geom.Box2I`

132 Current bounding box.

133 wcs : `lsst.afw.geom.SkyWcs

134 Current Wcs.

135

136 Returns

137 -------

138 finalBBox : `lsst.geom.Box2I`

139 Revised bounding box.

140 wcs : `lsst.afw.geom.SkyWcs`

141 Revised Wcs.

142 """

143 finalBBox = geom.Box2I(geom.Point2I(0, 0), bbox.getDimensions())

144 # shift the WCS by the same amount as the bbox; extra code is required

145 # because simply subtracting makes an Extent2I

146 pixPosOffset = geom.Extent2D(finalBBox.getMinX() - bbox.getMinX(),

147 finalBBox.getMinY() - bbox.getMinY())

148 wcs = wcs.copyAtShiftedPixelOrigin(pixPosOffset)

149 return finalBBox, wcs

150

151 def getSequentialPatchIndex(self, patchInfo):

152 """Return a single integer that uniquely identifies

153 the given patch within this tract.

154

155 Parameters

156 ----------

157 patchInfo : `lsst.skymap.PatchInfo`

158

159 Returns

160 -------

161 sequentialIndex : `int`

162 """

163 return self._tractBuilder.getSequentialPatchIndex(patchInfo)

164

165 def getSequentialPatchIndexFromPair(self, index):

166 """Return a single integer that uniquely identifies

167 the patch index within the tract.

168

169 Parameters

170 ----------

171 index : `lsst.skymap.Index2D`

172

173 Returns

174 -------

175 sequentialIndex : `int`

176 """

177 return self._tractBuilder.getSequentialPatchIndexFromPair(index)

178

179 def getPatchIndexPair(self, sequentialIndex):

180 """Convert sequential index into patch index (x,y) pair.

181

182 Parameters

183 ----------

184 sequentialIndex : `int`

185

186 Returns

187 -------

188 x, y : `lsst.skymap.Index2D`

189 """

190 return self._tractBuilder.getPatchIndexPair(sequentialIndex)

191

192 def findPatch(self, coord):

193 """Find the patch containing the specified coord.

194

195 Parameters

196 ----------

197 coord : `lsst.geom.SpherePoint`

198 ICRS sky coordinate to search for.

199

200 Returns

201 -------

202 result : `lsst.skymap.PatchInfo`

203 PatchInfo of patch whose inner bbox contains the specified coord

204

205 Raises

206 ------

207 LookupError

208 Raised if coord is not in tract or we cannot determine the

209 pixel coordinate (which likely means the coord is off the tract).

210 """

211 try:

212 pixel = self.wcs.skyToPixel(coord)

213 except (lsst.pex.exceptions.DomainError, lsst.pex.exceptions.RuntimeError):

214 # Point must be way off the tract

215 raise LookupError("Unable to determine pixel position for coordinate %s" % (coord,))

216 pixelInd = geom.Point2I(pixel)

217 if not self._bbox.contains(pixelInd):

218 raise LookupError("coord %s is not in tract %s" % (coord, self.tract_id))

219 # This should probably be the same as above because we only

220 # care about the INNER dimensions.

221 patchInd = tuple(int(pixelInd[i]/self.patch_inner_dimensions[i]) for i in range(2))

222 return self.getPatchInfo(patchInd)

223

224 def findPatchList(self, coordList):

225 """Find patches containing the specified list of coords.

226

227 Parameters

228 ----------

229 coordList : `list` of `lsst.geom.SpherePoint`

230 ICRS sky coordinates to search for.

231

232 Returns

233 -------

234 result : `list` of `lsst.skymap.PatchInfo`

235 List of PatchInfo for patches that contain, or may contain, the

236 specified region. The list will be empty if there is no overlap.

237

238 Notes

239 -----

240 **Warning:**

241

242 - This may give incorrect answers on regions that are larger than a

243 tract.

244

245 - This uses a naive algorithm that may find some patches that do not

246 overlap the region (especially if the region is not a rectangle

247 aligned along patch x,y).

248 """

249 box2D = geom.Box2D()

250 for coord in coordList:

251 try:

252 pixelPos = self.wcs.skyToPixel(coord)

253 except (lsst.pex.exceptions.DomainError, lsst.pex.exceptions.RuntimeError):

254 # The point is so far off the tract that its pixel position

255 # cannot be computed.

256 continue

257 box2D.include(pixelPos)

258 bbox = geom.Box2I(box2D)

259 bbox.grow(self.getPatchBorder())

260 bbox.clip(self._bbox)

261 if bbox.isEmpty():

262 return ()

263

264 llPatchInd = tuple(int(bbox.getMin()[i]/self.patch_inner_dimensions[i]) for i in range(2))

265 urPatchInd = tuple(int(bbox.getMax()[i]/self.patch_inner_dimensions[i]) for i in range(2))

266 return tuple(self.getPatchInfo((xInd, yInd))

267 for xInd in range(llPatchInd[0], urPatchInd[0]+1)

268 for yInd in range(llPatchInd[1], urPatchInd[1]+1))

269

270 def getBBox(self):

271 """Get bounding box of tract (as an geom.Box2I)

272 """

273 return geom.Box2I(self._bbox)

274

275 bbox = property(getBBox)

276

277 def getCtrCoord(self):

278 """Get ICRS sky coordinate of center of tract

279 (as an lsst.geom.SpherePoint)

280 """

281 return self._ctrCoord

282

283 ctr_coord = property(getCtrCoord)

284

285 def getId(self):

286 """Get ID of tract

287 """

288 return self._id

289

290 tract_id = property(getId)

291

292 def getNumPatches(self):

293 """Get the number of patches in x, y.

294

295 Returns

296 -------

297 result : `lsst.skymap.Index2D`

298 The number of patches in x, y

299 """

300 return self._numPatches

301

302 num_patches = property(getNumPatches)

303

304 def getPatchBorder(self):

305 return self._tractBuilder.getPatchBorder()

306

307 patch_border = property(getPatchBorder)

308

309 def getPatchInfo(self, index):

310 """Return information for the specified patch.

311

312 Parameters

313 ----------

314 index : `typing.NamedTuple` ['x': `int`, 'y': `int`]

315 Index of patch, as a pair of ints;

316 or a sequential index as returned by getSequentialPatchIndex;

317 negative values are not supported.

318

319 Returns

320 -------

321 result : `lsst.skymap.PatchInfo`

322 The patch info for that index.

323

324 Raises

325 ------

326 IndexError

327 Raised if index is out of range.

328 """

329 return self._tractBuilder.getPatchInfo(index, self._wcs)

330

331 def getPatchInnerDimensions(self):

332 """Get dimensions of inner region of the patches (all are the same)

333 """

334 return self._tractBuilder.getPatchInnerDimensions()

335

336 patch_inner_dimensions = property(getPatchInnerDimensions)

337

338 def getTractOverlap(self):

339 """Get minimum overlap of adjacent sky tracts.

340 """

341 return self._tractOverlap

342

343 tract_overlap = property(getTractOverlap)

344

345 def getVertexList(self):

346 """Get list of ICRS sky coordinates of vertices that define the

347 boundary of the inner region.

348

349 Notes

350 -----

351 **warning:** this is not a deep copy.

352 """

353 return self._vertexCoordList

354

355 vertex_list = property(getVertexList)

356

357 # TODO: Remove with DM-44799

358 @deprecated(reason="getInnerSkyPolygon()/inner_sky_polygon has been deprecated in favor of "

359 "inner_sky_region, and will be removed after v28.",

360 category=FutureWarning, version=28)

361 def getInnerSkyPolygon(self):

362 """Get inner on-sky region as a sphgeom.ConvexPolygon.

363 """

364 skyUnitVectors = [sp.getVector() for sp in self.getVertexList()]

365 return ConvexPolygon.convexHull(skyUnitVectors)

366

367 inner_sky_polygon = property(getInnerSkyPolygon)

368

369 def getInnerSkyRegion(self):

370 """Get inner on-sky region.

371 """

372 if self._innerBoxCorners:

373 return Box(point1=self._innerBoxCorners[0], point2=self._innerBoxCorners[1])

374 else:

375 skyUnitVectors = [sp.getVector() for sp in self.getVertexList()]

376 return ConvexPolygon.convexHull(skyUnitVectors)

377

378 inner_sky_region = property(getInnerSkyRegion)

379

380 def getOuterSkyPolygon(self):

381 """Get outer on-sky region as a sphgeom.ConvexPolygon

382 """

383 return makeSkyPolygonFromBBox(bbox=self.getBBox(), wcs=self.getWcs())

384

385 outer_sky_polygon = property(getOuterSkyPolygon)

386

387 def getWcs(self):

388 """Get WCS of tract.

389

390 Returns

391 -------

392 wcs : `lsst.afw.geom.SkyWcs`

393 The WCS of this tract

394 """

395 return self._wcs

396

397 wcs = property(getWcs)

398

399 def __str__(self):

400 return "TractInfo(id=%s)" % (self._id,)

401

402 def __repr__(self):

403 return "TractInfo(id=%s, ctrCoord=%s)" % (self._id, self._ctrCoord.getVector())

404

405 def __iter__(self):

406 xNum, yNum = self.getNumPatches()

407 for y in range(yNum):

408 for x in range(xNum):

409 yield self.getPatchInfo(Index2D(x=x, y=y))

410

411 def __len__(self):

412 xNum, yNum = self.getNumPatches()

413 return xNum*yNum

414

415 def __getitem__(self, index):

416 return self.getPatchInfo(index)

417

418 def contains(self, coord):

419 """Does this tract contain the coordinate?"""

420 try:

421 pixel = self.getWcs().skyToPixel(coord)

422 except (lsst.pex.exceptions.DomainError, lsst.pex.exceptions.RuntimeError):

423 # Point must be way off the tract

424 return False

425 if not np.isfinite(pixel.getX()) or not np.isfinite(pixel.getY()):

426 # Point is definitely off the tract

427 return False

428 return self.getBBox().contains(geom.Point2I(pixel))

429

430

431class ExplicitTractInfo(TractInfo):

432 """Information for a tract specified explicitly.

433

434 A tract is placed at the explicitly defined coordinates, with the nominated

435 radius. The tracts are square (i.e., the radius is really a half-size).

436

437 Parameters

438 ----------

439 id : : `int`

440 tract ID

441 tractBuilder : Subclass of `lsst.skymap.BaseTractBuilder`

442 Object used to compute patch geometry.

443 ctrCoord : `lsst.geom.SpherePoint`

444 ICRS sky coordinate of center of inner region of tract; also used as

445 the CRVAL for the WCS.

446 radius : `lsst.geom.Angle`

447 Radius of the tract.

448 tractOverlap : `lsst.geom.Angle`

449 Minimum overlap between adjacent sky tracts; this defines the minimum

450 distance the tract extends beyond the inner region in all directions.

451 wcs : `lsst.afw.image.SkyWcs`

452 WCS for tract. The reference pixel will be shifted as required so that

453 the lower left-hand pixel (index 0,0) has pixel position 0.0, 0.0.

454 innerBoxCorners : `list` [`lsst.sphgeom.LonLat`], optional

455 If set then the ``inner_sky_region`` will be a `lsst.sphgeom.Box` with

456 these corners as oppsed to a `lsst.sphgeom.ConvexPolygon` built from

457 the ``vertex_list``.

458 """

459 def __init__(self, id, tractBuilder, ctrCoord, radius, tractOverlap, wcs, innerBoxCorners=None):

460 # We don't want TractInfo setting the bbox on the basis of vertices,

461 # but on the radius.

462 vertexList = []

463 self._radius = radius

464 super(ExplicitTractInfo, self).__init__(

465 id,

466 tractBuilder,

467 ctrCoord,

468 vertexList,

469 tractOverlap,

470 wcs,

471 innerBoxCorners=innerBoxCorners,

472 )

473 # Shrink the box slightly to make sure the vertices are in the tract

474 bboxD = geom.BoxD(self.getBBox())

475 bboxD.grow(-0.001)

476 finalWcs = self.getWcs()

477 self._vertexCoordList = finalWcs.pixelToSky(bboxD.getCorners())

478

479 def _minimumBoundingBox(self, wcs):

480 """Calculate the minimum bounding box for the tract, given the WCS, and

481 the nominated radius.

482 """

483 bbox = geom.Box2D()

484 for i in range(4):

485 cornerCoord = self._ctrCoord.offset(i*90*geom.degrees, self._radius + self._tractOverlap)

486 pixPos = wcs.skyToPixel(cornerCoord)

487 bbox.include(pixPos)

488 return bbox

Coverage for python / lsst / skymap / tractInfo.py: 30%

150 statements