Coverage for python/lsst/skymap/baseSkyMap.py: 72%

1# LSST Data Management System

4# This product includes software developed by the

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

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

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

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

10# (at your option) any later version.

11#

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

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

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

15# GNU General Public License for more details.

16#

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

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

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

20#

22"""

23todo: Consider tweaking pixel scale so the average scale is as specified,

24rather than the scale at the center.

25"""

27__all__ = ["BaseSkyMapConfig", "BaseSkyMap"]

29import hashlib

30import numpy as np

32import lsst.geom as geom

33import lsst.pex.config as pexConfig

34from lsst.geom import Angle, arcseconds, degrees

35from . import detail

36from .tractBuilder import tractBuilderRegistry

39class BaseSkyMapConfig(pexConfig.Config):

40 tractBuilder = tractBuilderRegistry.makeField(

41 doc="Algorithm for creating patches within the tract.",

42 default="legacy"

43 )

45 tractOverlap = pexConfig.Field(

46 doc="minimum overlap between adjacent sky tracts, on the sky (deg)",

47 dtype=float,

48 default=1.0,

49 )

50 pixelScale = pexConfig.Field(

51 doc="nominal pixel scale (arcsec/pixel)",

52 dtype=float,

53 default=0.333

54 )

55 projection = pexConfig.Field(

56 doc="one of the FITS WCS projection codes, such as:"

57 "- STG: stereographic projection"

58 "- MOL: Molleweide's projection"

59 "- TAN: tangent-plane projection",

60 dtype=str,

61 default="STG",

62 )

63 rotation = pexConfig.Field(

64 doc="Rotation for WCS (deg)",

65 dtype=float,

66 default=0,

67 )

69 # Backwards compatibility

70 # We can't use the @property decorator because it makes pexConfig sad.

71 def getPatchInnerDimensions(self):

72 """Get the patch inner dimensions, for backwards compatibility.

74 This value is only used with the ``legacy`` tract builder,

75 and is ignored otherwise. In general, the config should be

76 accessed directly with

77 ``config.tractBuilder["legacy"].patchInnerDimensions``.

79 Returns

80 -------

81 innerDimensions : `list` [`int`, `int`]

82 """

83 return self.tractBuilder["legacy"].patchInnerDimensions

85 def setPatchInnerDimensions(self, value):

86 """Set the patch inner dimensions, for backwards compatibility.

88 This value is only used with the ``legacy`` tract builder,

89 and is ignored otherwise. In general, the config should be

90 accessed directly with

91 ``config.tractBuilder["legacy"].patchInnerDimensions``.

93 Parameters

94 ----------

95 value : `list` [`int`, `int`]

96 """

97 self.tractBuilder["legacy"].patchInnerDimensions = value

99 patchInnerDimensions = property(getPatchInnerDimensions, setPatchInnerDimensions)

100

101 def getPatchBorder(self):

102 """Get the patch border, for backwards compatibility.

103

104 This value is only used with the ``legacy`` tract builder,

105 and is ignored otherwise. In general, the config should be

106 accessed directly with config.tractBuilder["legacy"].patchBorder.

107

108 Returns

109 -------

110 border: `int`

111 """

112 return self.tractBuilder["legacy"].patchBorder

113

114 def setPatchBorder(self, value):

115 """Set the patch border, for backwards compatibility.

116

117 This value is only used with the ``legacy`` tract builder,

118 and is ignored otherwise. In general, the config should be

119 accessed directly with config.tractBuilder["legacy"].patchBorder.

120

121 Parameters

122 ----------

123 border: `int`

124 """

125 self.tractBuilder["legacy"].patchBorder = value

126

127 patchBorder = property(getPatchBorder, setPatchBorder)

128

129

130class BaseSkyMap:

131 """A collection of overlapping Tracts that map part or all of the sky.

132

133 See TractInfo for more information.

134

135 Parameters

136 ----------

137 config : `BaseSkyMapConfig` or `None` (optional)

138 The configuration for this SkyMap; if None use the default config.

139

140 Notes

141 -----

142 BaseSkyMap is an abstract base class. Subclasses must do the following:

143 define ``__init__`` and have it construct the TractInfo objects and put

144 them in ``__tractInfoList__`` define ``__getstate__`` and ``__setstate__``

145 to allow pickling (the butler saves sky maps using pickle);

146 see DodecaSkyMap for an example of how to do this. (Most of that code could

147 be moved into this base class, but that would make it harder to handle

148 older versions of pickle data.) define updateSha1 to add any

149 subclass-specific state to the hash.

150

151 All SkyMap subclasses must be conceptually immutable; they must always

152 refer to the same set of mathematical tracts and patches even if the in-

153 memory representation of those objects changes.

154 """

155

156 ConfigClass = BaseSkyMapConfig

157

158 def __init__(self, config=None):

159 if config is None:

160 config = self.ConfigClass()

161 config.freeze() # just to be sure, e.g. for pickling

162 self.config = config

163 self._tractInfoList = []

164 self._wcsFactory = detail.WcsFactory(

165 pixelScale=Angle(self.config.pixelScale, arcseconds),

166 projection=self.config.projection,

167 rotation=Angle(self.config.rotation, degrees),

168 )

169 self._sha1 = None

170 self._tractBuilder = config.tractBuilder.apply()

171

172 def findTract(self, coord):

173 """Find the tract whose center is nearest the specified coord.

174

175 Parameters

176 ----------

177 coord : `lsst.geom.SpherePoint`

178 ICRS sky coordinate to search for.

179

180 Returns

181 -------

182 result : `TractInfo`

183 TractInfo of tract whose center is nearest the specified coord.

184

185 Notes

186 -----

187 - If coord is equidistant between multiple sky tract centers then one

188 is arbitrarily chosen.

189

190 - The default implementation is not very efficient; subclasses may wish

191 to override.

192

193 .. warning::

194

195 If tracts do not cover the whole sky then the returned tract may not

196 include the coord.

197 """

198 distTractInfoList = []

199 for i, tractInfo in enumerate(self):

200 angSep = coord.separation(tractInfo.getCtrCoord()).asDegrees()

201 # include index in order to disambiguate identical angSep values

202 distTractInfoList.append((angSep, i, tractInfo))

203 distTractInfoList.sort()

204 return distTractInfoList[0][2]

205

206 def findTractIdArray(self, ra, dec, degrees=False):

207 """Find array of tract IDs with vectorized operations (where

208 supported).

209

210 If a given sky map does not support vectorized operations, then a loop

211 over findTract will be called.

212

213 Parameters

214 ----------

215 ra : `numpy.ndarray`

216 Array of Right Ascension. Units are radians unless

217 degrees=True.

218 dec : `numpy.ndarray`

219 Array of Declination. Units are radians unless

220 degrees=True.

221 degrees : `bool`, optional

222 Input ra, dec arrays are degrees if `True`.

223

224 Returns

225 -------

226 tractId : `numpy.ndarray`

227 Array of tract IDs

228

229 Notes

230 -----

231 - If coord is equidistant between multiple sky tract centers then one

232 is arbitrarily chosen.

233

234 .. warning::

235

236 If tracts do not cover the whole sky then the returned tract may not

237 include the given ra/dec.

238 """

239 units = geom.degrees if degrees else geom.radians

240 coords = [geom.SpherePoint(r*units, d*units) for r, d in zip(np.atleast_1d(ra),

241 np.atleast_1d(dec))]

242

243 tractId = np.array([self.findTract(coord).getId() for coord in coords])

244

245 return tractId

246

247 def findTractPatchList(self, coordList):

248 """Find tracts and patches that overlap a region.

249

250 Parameters

251 ----------

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

253 List of ICRS sky coordinates to search for.

254

255 Returns

256 -------

257 reList : `list` of (`TractInfo`, `list` of `PatchInfo`)

258 For tracts and patches that contain, or may contain, the specified

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

260

261 Notes

262 -----

263 .. warning::

264

265 This uses a naive algorithm that may find some tracts and patches

266 that do not overlap the region (especially if the region is not a

267 rectangle aligned along patch x, y).

268 """

269 retList = []

270 for tractInfo in self:

271 patchList = tractInfo.findPatchList(coordList)

272 if patchList:

273 retList.append((tractInfo, patchList))

274 return retList

275

276 def findClosestTractPatchList(self, coordList):

277 """Find closest tract and patches that overlap coordinates.

278

279 Parameters

280 ----------

281 coordList : `lsst.geom.SpherePoint`

282 List of ICRS sky coordinates to search for.

283

284 Returns

285 -------

286 retList : `list`

287 list of (TractInfo, list of PatchInfo) for tracts and patches

288 that contain, or may contain, the specified region.

289 The list will be empty if there is no overlap.

290 """

291 retList = []

292 for coord in coordList:

293 tractInfo = self.findTract(coord)

294 patchList = tractInfo.findPatchList(coordList)

295 if patchList and not (tractInfo, patchList) in retList: 295 ↛ 292line 295 didn't jump to line 292 because the condition on line 295 was always true

296 retList.append((tractInfo, patchList))

297 return retList

298

299 def findTractIdPatchIdArray(self, ra, dec, degrees=False):

300 """Find array of tract IDs and patch IDs with vectorized operations

301 (where supported).

302

303 If a given sky map does not support vectorized operations, then

304 loops will be called.

305

306 Parameters

307 ----------

308 ra : `numpy.ndarray`

309 Array of Right Ascension. Units are radians unless

310 degrees=True.

311 dec : `numpy.ndarray`

312 Array of Declination. Units are radians unless

313 degrees=True.

314 degrees : `bool`, optional

315 Input ra, dec arrays are degrees if `True`.

316

317 Returns

318 -------

319 tractId : `numpy.ndarray`

320 Array of tract IDs.

321 patchId : `numpy.ndarray`

322 Array of sequential patch IDs. -1 if there is no appropriate

323 patch.

324

325 Notes

326 -----

327 - If coord is equidistant between multiple sky tract centers then one

328 is arbitrarily chosen.

329

330 .. warning::

331

332 If tracts do not cover the whole sky then the returned tract may not

333 include the given ra/dec.

334 """

335 from scipy.ndimage import value_indices

336

337 _ra = np.atleast_1d(ra)

338 _dec = np.atleast_1d(dec)

339

340 # This will be vectorized if possible.

341 tractIds = self.findTractIdArray(_ra.ravel(), _dec.ravel(), degrees=degrees)

342

343 # This will be vectorized within each patch.

344 patchIds = np.zeros(len(tractIds), dtype=np.int32) - 1

345

346 inds = value_indices(tractIds)

347 for tractId, (indices,) in inds.items():

348 tract = self[tractId]

349 wcs = tract.wcs

350 x, y = wcs.skyToPixelArray(_ra[indices], _dec[indices], degrees=degrees)

351 xInd = np.floor(x).astype(np.int64) // tract.patch_inner_dimensions[0]

352 yInd = np.floor(y).astype(np.int64) // tract.patch_inner_dimensions[1]

353 nx, ny = tract.num_patches

354 patchIds[indices] = nx * yInd + xInd

355

356 # Check for overflows

357 bad = (xInd < 0) | (xInd >= nx) | (yInd < 0) | (yInd >= ny)

358 if bad.sum() > 0:

359 tractIds[indices[bad]] = -1

360 patchIds[indices[bad]] = -1

361

362 return (tractIds, patchIds)

363

364 def __getitem__(self, ind):

365 return self._tractInfoList[ind]

366

367 def __iter__(self):

368 return iter(self._tractInfoList)

369

370 def __len__(self):

371 return len(self._tractInfoList)

372

373 def __hash__(self):

374 return hash(self.getSha1())

375

376 def __eq__(self, other):

377 try:

378 return self.getSha1() == other.getSha1()

379 except AttributeError:

380 return NotImplemented

381

382 def __ne__(self, other):

383 return not (self == other)

384

385 def logSkyMapInfo(self, log):

386 """Write information about a sky map to supplied log

387

388 Parameters

389 ----------

390 log : `logging.Logger`

391 Log object that information about skymap will be written.

392 """

393 log.info("sky map has %s tracts" % (len(self),))

394 for tractInfo in self:

395 wcs = tractInfo.getWcs()

396 posBox = geom.Box2D(tractInfo.getBBox())

397 pixelPosList = (

398 posBox.getMin(),

399 geom.Point2D(posBox.getMaxX(), posBox.getMinY()),

400 posBox.getMax(),

401 geom.Point2D(posBox.getMinX(), posBox.getMaxY()),

402 )

403 skyPosList = [wcs.pixelToSky(pos).getPosition(geom.degrees) for pos in pixelPosList]

404 posStrList = ["(%0.3f, %0.3f)" % tuple(skyPos) for skyPos in skyPosList]

405 log.info("tract %s has corners %s (RA, Dec deg) and %s x %s patches" %

406 (tractInfo.getId(), ", ".join(posStrList),

407 tractInfo.getNumPatches()[0], tractInfo.getNumPatches()[1]))

408

409 def getSha1(self):

410 """Return a SHA1 hash that uniquely identifies this SkyMap instance.

411

412 Returns

413 -------

414 sha1 : `bytes`

415 A 20-byte hash that uniquely identifies this SkyMap instance.

416

417 Notes

418 -----

419 Subclasses should almost always override ``updateSha1`` instead of

420 this function to add subclass-specific state to the hash.

421 """

422 if self._sha1 is None:

423 sha1 = hashlib.sha1()

424 sha1.update(type(self).__name__.encode('utf-8'))

425 configPacked = self._tractBuilder.getPackedConfig(self.config)

426 sha1.update(configPacked)

427 self.updateSha1(sha1)

428 self._sha1 = sha1.digest()

429 return self._sha1

430

431 def updateSha1(self, sha1):

432 """Add subclass-specific state or configuration options to the SHA1.

433

434 Parameters

435 ----------

436 sha1 : `hashlib.sha1`

437 A hashlib object on which `update` can be called to add

438 additional state to the hash.

439

440 Notes

441 -----

442 This method is conceptually "protected" : it should be reimplemented by

443 all subclasses, but called only by the base class implementation of

444 `getSha1` .

445 """

446 raise NotImplementedError()

447

448 SKYMAP_RUN_COLLECTION_NAME = "skymaps"

449

450 SKYMAP_DATASET_TYPE_NAME = "skyMap"

451

452 def register(self, name, butler):

453 """Add skymap, tract, and patch Dimension entries to the given Gen3

454 Butler.

455

456 Parameters

457 ----------

458 name : `str`

459 The name of the skymap.

460 butler : `lsst.daf.butler.Butler`

461 The butler to add to.

462

463 Raises

464 ------

465 lsst.daf.butler.registry.ConflictingDefinitionError

466 Raised if a different skymap exists with the same name.

467

468 Notes

469 -----

470 Registering the same skymap multiple times (with the exact same

471 definition) is safe, but inefficient; most of the work of computing

472 the rows to be inserted must be done first in order to check for

473 consistency between the new skymap and any existing one.

474

475 Re-registering a skymap with different tract and/or patch definitions

476 but the same summary information may not be detected as a conflict but

477 will never result in updating the skymap; there is intentionally no

478 way to modify a registered skymap (aside from manual administrative

479 operations on the database), as it is hard to guarantee that this can

480 be done without affecting reproducibility.

481 """

482 numPatches = [tractInfo.getNumPatches() for tractInfo in self]

483 nxMax = max(nn[0] for nn in numPatches)

484 nyMax = max(nn[1] for nn in numPatches)

485

486 skyMapRecord = {

487 "skymap": name,

488 "hash": self.getSha1(),

489 "tract_max": len(self),

490 "patch_nx_max": nxMax,

491 "patch_ny_max": nyMax,

492 }

493 butler.registry.registerRun(self.SKYMAP_RUN_COLLECTION_NAME)

494 # Kind of crazy that we've got three different capitalizations of

495 # "skymap" here, but that's what the various conventions (or at least

496 # precedents) dictate.

497 from lsst.daf.butler import DatasetType

498 from lsst.daf.butler.registry import ConflictingDefinitionError

499 datasetType = DatasetType(

500 name=self.SKYMAP_DATASET_TYPE_NAME,

501 dimensions=["skymap"],

502 storageClass="SkyMap",

503 universe=butler.dimensions

504 )

505 butler.registry.registerDatasetType(datasetType)

506 with butler.transaction():

507 try:

508 inserted = butler.registry.syncDimensionData("skymap", skyMapRecord)

509 except ConflictingDefinitionError as err:

510 raise ConflictingDefinitionError(

511 f"SkyMap with hash {self.getSha1().hex()} is already registered with a different name."

512 ) from err

513 if inserted:

514 for tractInfo in self:

515 tractId = tractInfo.getId()

516 tractRegion = tractInfo.getOuterSkyPolygon()

517 tractWcs = tractInfo.getWcs()

518 tractRecord = dict(

519 skymap=name,

520 tract=tractId,

521 region=tractRegion,

522 )

523 butler.registry.insertDimensionData("tract", tractRecord)

524

525 patchRecords = []

526 for patchInfo in tractInfo:

527 xx, yy = patchInfo.getIndex()

528 patchRecords.append(

529 dict(

530 skymap=name,

531 tract=tractId,

532 patch=tractInfo.getSequentialPatchIndex(patchInfo),

533 cell_x=xx,

534 cell_y=yy,

535 region=patchInfo.getOuterSkyPolygon(tractWcs),

536 )

537 )

538 butler.registry.insertDimensionData("patch", *patchRecords)

539

540 butler.put(self, datasetType, {"skymap": name}, run=self.SKYMAP_RUN_COLLECTION_NAME)