Coverage for tests / test_gbdesAstrometricFit.py: 9%
604 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-24 08:49 +0000
1# This file is part of drp_tasks
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import os.path
23import unittest
24from copy import copy
26import astropy.table
27import astropy.time
28import astropy.units as u
29import numpy as np
30import pandas as pd
31import wcsfit
32import yaml
33from astropy.coordinates import Distance, SkyCoord
34from smatch.matcher import Matcher
36import lsst.afw.geom as afwgeom
37import lsst.afw.table as afwTable
38import lsst.geom
39import lsst.utils
40from lsst import sphgeom
41from lsst.daf.base import PropertyList
42from lsst.drp.tasks.gbdesAstrometricFit import (
43 GbdesAstrometricFitConfig,
44 GbdesAstrometricFitTask,
45 GbdesGlobalAstrometricFitConfig,
46 GbdesGlobalAstrometricFitTask,
47 calculate_apparent_motion,
48)
49from lsst.meas.algorithms import ReferenceObjectLoader
50from lsst.meas.algorithms.testUtils import MockRefcatDataId
51from lsst.pipe.base import InMemoryDatasetHandle
53TESTDIR = os.path.abspath(os.path.dirname(__file__))
56class TestGbdesAstrometricFit(lsst.utils.tests.TestCase):
57 """This class tests `GbdesAstrometricFit` using real `visitSummaryTable`s
58 from HSC RC2 processing, with simulated sources for those visits.
59 """
61 @classmethod
62 def setUpClass(cls):
63 # Set random seed
64 np.random.seed(1234)
66 # Make fake data
67 cls.datadir = os.path.join(TESTDIR, "data")
69 cls.fieldNumber = 0
70 cls.nFields = 1
71 cls.instrumentName = "HSC"
72 cls.refEpoch = 57205.5
74 # Make test inputVisitSummary. VisitSummaryTables are taken from
75 # collection HSC/runs/RC2/w_2022_20/DM-34794
76 cls.testVisits = [1176, 17900, 17930, 17934]
77 cls.inputVisitSummary = []
78 for testVisit in cls.testVisits:
79 visSum = afwTable.ExposureCatalog.readFits(
80 os.path.join(cls.datadir, f"visitSummary_{testVisit}.fits")
81 )
82 cls.inputVisitSummary.append(visSum)
84 cls.config = GbdesAstrometricFitConfig()
85 cls.config.systematicError = 0
86 cls.config.devicePolyOrder = 4
87 cls.config.exposurePolyOrder = 6
88 cls.config.fitReserveFraction = 0
89 cls.config.fitReserveRandomSeed = 1234
90 cls.config.saveModelParams = True
91 cls.config.allowSelfMatches = True
92 cls.config.saveCameraModel = True
93 cls.config.applyRefCatProperMotion = False
94 cls.task = GbdesAstrometricFitTask(config=cls.config)
96 (
97 cls.exposureInfo,
98 cls.exposuresHelper,
99 cls.extensionInfo,
100 cls.instruments,
101 ) = cls.task._get_exposure_info(cls.inputVisitSummary, refEpoch=cls.refEpoch)
103 cls.fields, cls.fieldCenter, cls.fieldRadius = cls.task._prep_sky(
104 cls.inputVisitSummary, cls.exposureInfo.medianEpoch
105 )
107 # Bounding box of observations:
108 raMins, raMaxs = [], []
109 decMins, decMaxs = [], []
110 for visSum in cls.inputVisitSummary:
111 raMins.append(visSum["raCorners"].min())
112 raMaxs.append(visSum["raCorners"].max())
113 decMins.append(visSum["decCorners"].min())
114 decMaxs.append(visSum["decCorners"].max())
115 raMin = min(raMins)
116 raMax = max(raMaxs)
117 decMin = min(decMins)
118 decMax = max(decMaxs)
120 corners = [
121 lsst.geom.SpherePoint(raMin, decMin, lsst.geom.degrees).getVector(),
122 lsst.geom.SpherePoint(raMax, decMin, lsst.geom.degrees).getVector(),
123 lsst.geom.SpherePoint(raMax, decMax, lsst.geom.degrees).getVector(),
124 lsst.geom.SpherePoint(raMin, decMax, lsst.geom.degrees).getVector(),
125 ]
126 cls.boundingPolygon = sphgeom.ConvexPolygon(corners)
128 # Make random set of data in a bounding box determined by input visits
129 # Make wcs objects for the "true" model
130 cls.nStars, starIds, starRAs, starDecs = cls._make_simulated_stars(
131 raMin, raMax, decMin, decMax, cls.config.matchRadius
132 )
134 # Fraction of simulated stars in the reference catalog and science
135 # exposures
136 inReferenceFraction = 1
137 inScienceFraction = 1
139 # Make a reference catalog and load it into ReferenceObjectLoader
140 refDataId, deferredRefCat = cls._make_refCat(
141 starIds, starRAs, starDecs, inReferenceFraction, cls.boundingPolygon
142 )
143 cls.refObjectLoader = ReferenceObjectLoader([refDataId], [deferredRefCat])
144 cls.refObjectLoader.config.requireProperMotion = False
145 cls.refObjectLoader.config.anyFilterMapsToThis = "test_filter"
147 cls.task.refObjectLoader = cls.refObjectLoader
149 # Get True WCS for stars:
150 with open(os.path.join(cls.datadir, "sample_wcs.yaml"), "r") as f:
151 cls.trueModel = yaml.load(f, Loader=yaml.Loader)
153 trueWCSs = cls._make_wcs(cls.trueModel, cls.inputVisitSummary)
155 # Make source catalogs:
156 cls.inputCatalogRefs = cls._make_sourceCat(starIds, starRAs, starDecs, trueWCSs, inScienceFraction)
158 cls.colorCatalog = cls._make_colors(starRAs, starDecs)
160 cls.outputs = cls.task.run(
161 cls.inputCatalogRefs,
162 cls.inputVisitSummary,
163 instrumentName=cls.instrumentName,
164 refEpoch=cls.refEpoch,
165 refObjectLoader=cls.refObjectLoader,
166 )
168 @staticmethod
169 def _make_simulated_stars(raMin, raMax, decMin, decMax, matchRadius, nStars=10000):
170 """Generate random positions for "stars" in an RA/Dec box.
172 Parameters
173 ----------
174 raMin : `float`
175 Minimum RA for simulated stars.
176 raMax : `float`
177 Maximum RA for simulated stars.
178 decMin : `float`
179 Minimum Dec for simulated stars.
180 decMax : `float`
181 Maximum Dec for simulated stars.
182 matchRadius : `float`
183 Minimum allowed distance in arcsec between stars.
184 nStars : `int`, optional
185 Number of stars to simulate. Final number will be lower if any
186 too-close stars are dropped.
188 Returns
189 -------
190 nStars : `int`
191 Number of stars simulated.
192 starIds: `np.ndarray`
193 Unique identification number for stars.
194 starRAs: `np.ndarray`
195 Simulated Right Ascensions.
196 starDecs: `np.ndarray`
197 Simulated Declination.
198 """
199 starIds = np.arange(nStars)
200 starRAs = np.random.random(nStars) * (raMax - raMin) + raMin
201 starDecs = np.random.random(nStars) * (decMax - decMin) + decMin
202 # Remove neighbors:
203 with Matcher(starRAs, starDecs) as matcher:
204 idx = matcher.query_groups(matchRadius / 3600.0, min_match=2)
205 if len(idx) > 0:
206 neighbors = np.unique(np.concatenate(idx))
207 starRAs = np.delete(starRAs, neighbors)
208 starDecs = np.delete(starDecs, neighbors)
209 nStars = len(starRAs)
210 starIds = np.arange(nStars)
212 return nStars, starIds, starRAs, starDecs
214 @classmethod
215 def _make_refCat(cls, starIds, starRas, starDecs, inReferenceFraction, bounds):
216 """Make reference catalog from a subset of the simulated data
218 Parameters
219 ----------
220 starIds : `np.ndarray` [`int`]
221 Source ids for the simulated stars
222 starRas : `np.ndarray` [`float`]
223 RAs of the simulated stars
224 starDecs : `np.ndarray` [`float`]
225 Decs of the simulated stars
226 inReferenceFraction : float
227 Percentage of simulated stars to include in reference catalog
228 bounds : `lsst.sphgeom.ConvexPolygon`
229 Boundary of the reference catalog region
231 Returns
232 -------
233 refDataId : `lsst.meas.algorithms.testUtils.MockRefcatDataId`
234 Object that replicates the functionality of a dataId.
235 deferredRefCat : `lsst.pipe.base.InMemoryDatasetHandle`
236 Dataset handle for reference catalog.
237 """
238 nRefs = int(cls.nStars * inReferenceFraction)
239 refStarIndices = np.random.choice(cls.nStars, nRefs, replace=False)
240 # Make simpleCatalog to hold data, create datasetRef with `region`
241 # determined by bounding box used in above simulate.
242 refSchema = afwTable.SimpleTable.makeMinimalSchema()
243 idKey = refSchema.addField("sourceId", type="I")
244 fluxKey = refSchema.addField("test_filter_flux", units="nJy", type=np.float64)
245 raErrKey = refSchema.addField("coord_raErr", units="rad", type=np.float64)
246 decErrKey = refSchema.addField("coord_decErr", units="rad", type=np.float64)
247 pmraErrKey = refSchema.addField("pm_raErr", units="rad2 / yr", type=np.float64)
248 pmdecErrKey = refSchema.addField("pm_decErr", units="rad2 / yr", type=np.float64)
249 refCat = afwTable.SimpleCatalog(refSchema)
250 ref_md = PropertyList()
251 ref_md.set("REFCAT_FORMAT_VERSION", 1)
252 refCat.table.setMetadata(ref_md)
253 for i in refStarIndices:
254 record = refCat.addNew()
255 record.set(idKey, starIds[i])
256 record.setRa(lsst.geom.Angle(starRas[i], lsst.geom.degrees))
257 record.setDec(lsst.geom.Angle(starDecs[i], lsst.geom.degrees))
258 record.set(fluxKey, 1)
259 record.set(raErrKey, 0.00001)
260 record.set(decErrKey, 0.00001)
261 record.set(pmraErrKey, 1e-9)
262 record.set(pmdecErrKey, 1e-9)
263 refDataId = MockRefcatDataId(bounds)
264 deferredRefCat = InMemoryDatasetHandle(refCat, storageClass="SourceCatalog", htm7="mockRefCat")
266 return refDataId, deferredRefCat
268 @classmethod
269 def _make_sourceCat(cls, starIds, starRas, starDecs, trueWCSs, inScienceFraction):
270 """Make a `pd.DataFrame` catalog with the columns needed for the
271 object selector.
273 Parameters
274 ----------
275 starIds : `np.ndarray` [`int`]
276 Source ids for the simulated stars
277 starRas : `np.ndarray` [`float`]
278 RAs of the simulated stars
279 starDecs : `np.ndarray` [`float`]
280 Decs of the simulated stars
281 trueWCSs : `list` [`lsst.afw.geom.SkyWcs`]
282 WCS with which to simulate the source pixel coordinates
283 inReferenceFraction : float
284 Percentage of simulated stars to include in reference catalog
286 Returns
287 -------
288 sourceCat : `list` [`lsst.pipe.base.InMemoryDatasetHandle`]
289 List of reference to source catalogs.
290 """
291 inputCatalogRefs = []
292 # Take a subset of the simulated data
293 # Use true wcs objects to put simulated data into ccds
294 bbox = lsst.geom.BoxD(
295 lsst.geom.Point2D(
296 cls.inputVisitSummary[0][0]["bbox_min_x"], cls.inputVisitSummary[0][0]["bbox_min_y"]
297 ),
298 lsst.geom.Point2D(
299 cls.inputVisitSummary[0][0]["bbox_max_x"], cls.inputVisitSummary[0][0]["bbox_max_y"]
300 ),
301 )
302 bboxCorners = bbox.getCorners()
303 cls.inputCatalogRefs = []
304 for v, visit in enumerate(cls.testVisits):
305 nVisStars = int(cls.nStars * inScienceFraction)
306 visitStarIndices = np.random.choice(cls.nStars, nVisStars, replace=False)
307 visitStarIds = starIds[visitStarIndices]
308 visitStarRas = starRas[visitStarIndices]
309 visitStarDecs = starDecs[visitStarIndices]
310 sourceCats = []
311 for detector in trueWCSs[visit]:
312 detWcs = detector.getWcs()
313 detectorId = detector["id"]
314 radecCorners = detWcs.pixelToSky(bboxCorners)
315 detectorFootprint = sphgeom.ConvexPolygon([rd.getVector() for rd in radecCorners])
316 detectorIndices = detectorFootprint.contains(
317 (visitStarRas * u.degree).to(u.radian), (visitStarDecs * u.degree).to(u.radian)
318 )
319 nDetectorStars = detectorIndices.sum()
320 detectorArray = np.ones(nDetectorStars, dtype=bool) * detector["id"]
322 ones_like = np.ones(nDetectorStars)
323 zeros_like = np.zeros(nDetectorStars, dtype=bool)
325 x, y = detWcs.skyToPixelArray(
326 visitStarRas[detectorIndices], visitStarDecs[detectorIndices], degrees=True
327 )
329 origWcs = (cls.inputVisitSummary[v][cls.inputVisitSummary[v]["id"] == detectorId])[0].getWcs()
330 inputRa, inputDec = origWcs.pixelToSkyArray(x, y, degrees=True)
332 sourceDict = {}
333 sourceDict["detector"] = detectorArray
334 sourceDict["sourceId"] = visitStarIds[detectorIndices]
335 sourceDict["x"] = x
336 sourceDict["y"] = y
337 sourceDict["xErr"] = 1e-3 * ones_like
338 sourceDict["yErr"] = 1e-3 * ones_like
339 sourceDict["inputRA"] = inputRa
340 sourceDict["inputDec"] = inputDec
341 sourceDict["trueRA"] = visitStarRas[detectorIndices]
342 sourceDict["trueDec"] = visitStarDecs[detectorIndices]
343 for key in ["apFlux_12_0_flux", "apFlux_12_0_instFlux", "ixx", "iyy"]:
344 sourceDict[key] = ones_like
345 for key in [
346 "pixelFlags_edge",
347 "pixelFlags_saturated",
348 "pixelFlags_interpolatedCenter",
349 "pixelFlags_interpolated",
350 "pixelFlags_crCenter",
351 "pixelFlags_bad",
352 "hsmPsfMoments_flag",
353 "apFlux_12_0_flag",
354 "extendedness",
355 "sizeExtendedness",
356 "parentSourceId",
357 "deblend_nChild",
358 "ixy",
359 ]:
360 sourceDict[key] = zeros_like
361 sourceDict["apFlux_12_0_instFluxErr"] = 1e-2 * ones_like
362 sourceDict["detect_isPrimary"] = ones_like.astype(bool)
364 sourceCat = pd.DataFrame(sourceDict)
365 sourceCats.append(sourceCat)
367 visitSourceTable = pd.concat(sourceCats)
369 inputCatalogRef = InMemoryDatasetHandle(
370 visitSourceTable, storageClass="DataFrame", dataId={"visit": visit}
371 )
373 inputCatalogRefs.append(inputCatalogRef)
375 return inputCatalogRefs
377 @classmethod
378 def _make_wcs(cls, model, inputVisitSummaries):
379 """Make a `lsst.afw.geom.SkyWcs` from given model parameters
381 Parameters
382 ----------
383 model : `dict`
384 Dictionary with WCS model parameters
385 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
386 Visit summary catalogs
387 Returns
388 -------
389 catalogs : `dict` [`int`, `lsst.afw.table.ExposureCatalog`]
390 Visit summary catalogs with WCS set to input model
391 """
393 # Pixels will need to be rescaled before going into the mappings
394 xscale = inputVisitSummaries[0][0]["bbox_max_x"] - inputVisitSummaries[0][0]["bbox_min_x"]
395 yscale = inputVisitSummaries[0][0]["bbox_max_y"] - inputVisitSummaries[0][0]["bbox_min_y"]
397 catalogs = {}
398 schema = lsst.afw.table.ExposureTable.makeMinimalSchema()
399 schema.addField("visit", type="L", doc="Visit number")
400 for visitSum in inputVisitSummaries:
401 visit = visitSum[0]["visit"]
402 visitMapName = f"{visit}/poly"
403 visitModel = model[visitMapName]
405 catalog = lsst.afw.table.ExposureCatalog(schema)
406 catalog.resize(len(visitSum))
407 catalog["visit"] = visit
409 raDec = visitSum[0].getVisitInfo().getBoresightRaDec()
411 visitMapType = visitModel["Type"]
412 visitDict = {"Type": visitMapType}
413 if visitMapType == "Poly":
414 mapCoefficients = visitModel["XPoly"]["Coefficients"] + visitModel["YPoly"]["Coefficients"]
415 visitDict["Coefficients"] = mapCoefficients
417 for d, detector in enumerate(visitSum):
418 detectorId = detector["id"]
419 filterName = detector["physical_filter"]
420 detectorMapName = f"{filterName}/{detectorId}/poly"
421 detectorModel = model[detectorMapName]
423 detectorMapType = detectorModel["Type"]
424 mapDict = {detectorMapName: {"Type": detectorMapType}, visitMapName: visitDict}
425 if detectorMapType == "Poly":
426 mapCoefficients = (
427 detectorModel["XPoly"]["Coefficients"] + detectorModel["YPoly"]["Coefficients"]
428 )
429 mapDict[detectorMapName]["Coefficients"] = mapCoefficients
431 outWCS = cls.task._make_afw_wcs(
432 mapDict,
433 raDec.getRa(),
434 raDec.getDec(),
435 doNormalizePixels=True,
436 xScale=xscale,
437 yScale=yscale,
438 )
439 catalog[d].setId(detectorId)
440 catalog[d].setWcs(outWCS)
442 catalog.sort()
443 catalogs[visit] = catalog
445 return catalogs
447 @classmethod
448 def _make_colors(cls, starRas, starDecs):
449 """Make a catalog with the star magnitudes.
451 Parameters
452 ----------
453 starRas : `np.ndarray` [`float`]
454 RAs of the simulated stars
455 starDecs : `np.ndarray` [`float`]
456 Decs of the simulated stars
458 Returns
459 -------
460 colorCatalog : `lsst.afw.table.SimpleCatalog`
461 Catalog with star magnitudes.
462 """
463 bands = ["g", "r", "i", "z", "y"]
464 nStars = len(starRas)
466 # Make a catalog following what is done in `fgcmCal`.
467 schema = afwTable.SimpleTable.makeMinimalSchema()
468 schema.addField(
469 "mag_std_noabs",
470 type="ArrayF",
471 doc="Standard magnitude (no absolute calibration)",
472 size=len(bands),
473 )
474 colorCatalog = afwTable.SimpleCatalog(schema)
475 colorCatalog.resize(len(starRas))
476 colorCatalog["coord_ra"] = (starRas * u.degree + 10 * np.random.randn(nStars) * u.mas).to(u.radian)
477 colorCatalog["coord_dec"] = (starDecs * u.degree + 10 * np.random.randn(nStars) * u.mas).to(u.radian)
479 magMin = 19
480 magMax = 23
481 for i in range(len(bands)):
482 colorCatalog["mag_std_noabs"][:, i] = np.random.random(nStars) * (magMax - magMin) + magMin
484 md = PropertyList()
485 md.set("BANDS", bands)
486 colorCatalog.setMetadata(md)
487 return colorCatalog
489 def test_get_exposure_info(self):
490 """Test that information for input exposures is as expected and that
491 the WCS in the class object gives approximately the same results as the
492 input `lsst.afw.geom.SkyWcs`.
493 """
495 # The total number of extensions is the number of detectors for each
496 # visit plus one for the reference catalog for each field.
497 totalExtensions = sum([len(visSum) for visSum in self.inputVisitSummary]) + self.nFields
499 self.assertEqual(totalExtensions, len(self.extensionInfo.visit))
501 taskVisits = set(self.extensionInfo.visit)
502 refVisits = np.arange(0, -1 * self.nFields, -1).tolist()
503 self.assertEqual(taskVisits, set(self.testVisits + refVisits))
505 xx = np.linspace(0, 2000, 3)
506 yy = np.linspace(0, 4000, 6)
507 xgrid, ygrid = np.meshgrid(xx, yy)
508 for visSum in self.inputVisitSummary:
509 visit = visSum[0]["visit"]
510 for detectorInfo in visSum:
511 detector = detectorInfo["id"]
512 extensionIndex = np.flatnonzero(
513 (self.extensionInfo.visit == visit) & (self.extensionInfo.detector == detector)
514 )[0]
515 fitWcs = self.extensionInfo.wcs[extensionIndex]
516 calexpWcs = detectorInfo.getWcs()
518 tanPlaneXY = np.array([fitWcs.toWorld(x, y) for (x, y) in zip(xgrid.ravel(), ygrid.ravel())])
520 calexpra, calexpdec = calexpWcs.pixelToSkyArray(xgrid.ravel(), ygrid.ravel(), degrees=True)
522 # The pixel origin maps to a position slightly off from the
523 # tangent plane origin, so we want to use this value for the
524 # tangent-plane-to-sky part of the mapping.
525 tangentPlaneCenter = fitWcs.toWorld(
526 calexpWcs.getPixelOrigin().getX(), calexpWcs.getPixelOrigin().getY()
527 )
528 tangentPlaneOrigin = lsst.geom.Point2D(tangentPlaneCenter)
529 skyOrigin = calexpWcs.pixelToSky(
530 calexpWcs.getPixelOrigin().getX(), calexpWcs.getPixelOrigin().getY()
531 )
532 cdMatrix = afwgeom.makeCdMatrix(1.0 * lsst.geom.degrees, 0.0 * lsst.geom.degrees, True)
533 iwcToSkyWcs = afwgeom.makeSkyWcs(tangentPlaneOrigin, skyOrigin, cdMatrix)
534 newRAdeg, newDecdeg = iwcToSkyWcs.pixelToSkyArray(
535 tanPlaneXY[:, 0], tanPlaneXY[:, 1], degrees=True
536 )
538 np.testing.assert_allclose(calexpra, newRAdeg)
539 np.testing.assert_allclose(calexpdec, newDecdeg)
541 def test_refCatLoader(self):
542 """Test that we can load objects from refCat"""
544 tmpAssociations = wcsfit.FoFClass(
545 self.fields,
546 self.instruments,
547 self.exposuresHelper,
548 [self.fieldRadius.asDegrees()],
549 (self.task.config.matchRadius * u.arcsec).to(u.degree).value,
550 )
552 self.task._load_refcat(
553 self.refObjectLoader,
554 self.extensionInfo,
555 associations=tmpAssociations,
556 center=self.fieldCenter,
557 radius=self.fieldRadius,
558 epoch=self.exposureInfo.medianEpoch,
559 )
561 # We have only loaded one catalog, so getting the 'matches' should just
562 # return the same objects we put in, except some random objects that
563 # are too close together.
564 tmpAssociations.sortMatches(self.fieldNumber, minMatches=1)
566 nMatches = (np.array(tmpAssociations.sequence) == 0).sum()
568 self.assertLessEqual(nMatches, self.nStars)
569 self.assertGreater(nMatches, self.nStars * 0.9)
571 def test_loading_and_association(self):
572 """Test that objects can be loaded and correctly associated."""
573 # Running `_load_catalogs_and_associate` changes the input WCSs, so
574 # recalculate them here so that the variables shared among tests are
575 # not affected.
576 _, exposuresHelper, extensionInfo, instruments = self.task._get_exposure_info(
577 self.inputVisitSummary, refEpoch=self.refEpoch
578 )
580 tmpAssociations = wcsfit.FoFClass(
581 self.fields,
582 self.instruments,
583 exposuresHelper,
584 [self.fieldRadius.asDegrees()],
585 (self.task.config.matchRadius * u.arcsec).to(u.degree).value,
586 )
587 self.task._load_catalogs_and_associate(tmpAssociations, self.inputCatalogRefs, extensionInfo)
589 tmpAssociations.sortMatches(self.fieldNumber, minMatches=2)
591 matchIds = []
592 correctMatches = []
593 for s, e, o in zip(tmpAssociations.sequence, tmpAssociations.extn, tmpAssociations.obj):
594 objVisitInd = extensionInfo.visitIndex[e]
595 objDet = extensionInfo.detector[e]
596 extnInds = self.inputCatalogRefs[objVisitInd].get()["detector"] == objDet
597 objInfo = self.inputCatalogRefs[objVisitInd].get()[extnInds].iloc[o]
598 if s == 0:
599 if len(matchIds) > 0:
600 correctMatches.append(len(set(matchIds)) == 1)
601 matchIds = []
603 matchIds.append(objInfo["sourceId"])
605 # A few matches may incorrectly associate sources because of the random
606 # positions
607 self.assertGreater(sum(correctMatches), len(correctMatches) * 0.95)
609 def test_make_outputs(self):
610 """Test that the run method recovers the input model parameters."""
611 for v, visit in enumerate(self.testVisits):
612 visitSummary = self.inputVisitSummary[v]
613 outputWcsCatalog = self.outputs.outputWcss[visit]
614 visitSources = self.inputCatalogRefs[v].get()
615 for d, detectorRow in enumerate(visitSummary):
616 detectorId = detectorRow["id"]
617 fitwcs = outputWcsCatalog[d].getWcs()
618 detSources = visitSources[visitSources["detector"] == detectorId]
619 fitRA, fitDec = fitwcs.pixelToSkyArray(detSources["x"], detSources["y"], degrees=True)
620 dRA = fitRA - detSources["trueRA"]
621 dDec = fitDec - detSources["trueDec"]
622 # Check that input coordinates match the output coordinates
623 self.assertAlmostEqual(np.mean(dRA), 0)
624 self.assertAlmostEqual(np.std(dRA), 0)
625 self.assertAlmostEqual(np.mean(dDec), 0)
626 self.assertAlmostEqual(np.std(dDec), 0)
628 def test_compute_model_params(self):
629 """Test the optional model parameters and covariance output."""
630 modelParams = pd.DataFrame(self.outputs.modelParams)
631 # Check that DataFrame is the expected size.
632 shape = modelParams.shape
633 self.assertEqual(shape[0] + 4, shape[1])
634 # Check that covariance matrix is symmetric.
635 covariance = (modelParams.iloc[:, 4:]).to_numpy()
636 np.testing.assert_allclose(covariance, covariance.T, atol=1e-18)
638 def test_run(self):
639 """Test that run method recovers the input model parameters"""
640 outputMaps = self.outputs.fitModel.mapCollection.getParamDict()
642 for v, visit in enumerate(self.testVisits):
643 visitSummary = self.inputVisitSummary[v]
644 visitMapName = f"{visit}/poly"
646 origModel = self.trueModel[visitMapName]
647 if origModel["Type"] != "Identity":
648 fitModel = outputMaps[visitMapName]
649 origXPoly = origModel["XPoly"]["Coefficients"]
650 origYPoly = origModel["YPoly"]["Coefficients"]
651 fitXPoly = fitModel[: len(origXPoly)]
652 fitYPoly = fitModel[len(origXPoly) :]
654 absDiffX = abs(fitXPoly - origXPoly)
655 absDiffY = abs(fitYPoly - origYPoly)
656 # Check that input visit model matches fit
657 np.testing.assert_array_less(absDiffX, 1e-6)
658 np.testing.assert_array_less(absDiffY, 1e-6)
659 for d, detectorRow in enumerate(visitSummary):
660 detectorId = detectorRow["id"]
661 filterName = detectorRow["physical_filter"]
662 detectorMapName = f"{filterName}/{detectorId}/poly"
663 origModel = self.trueModel[detectorMapName]
664 if (origModel["Type"] != "Identity") and (v == 0):
665 fitModel = outputMaps[detectorMapName]
666 origXPoly = origModel["XPoly"]["Coefficients"]
667 origYPoly = origModel["YPoly"]["Coefficients"]
668 fitXPoly = fitModel[: len(origXPoly)]
669 fitYPoly = fitModel[len(origXPoly) :]
670 absDiffX = abs(fitXPoly - origXPoly)
671 absDiffY = abs(fitYPoly - origYPoly)
672 # Check that input detector model matches fit
673 np.testing.assert_array_less(absDiffX, 1e-7)
674 np.testing.assert_array_less(absDiffY, 1e-7)
676 def test_missingWcs(self):
677 """Test that task does not fail when the input WCS is None for one
678 extension and that the fit WCS for that extension returns a finite
679 result.
680 """
681 # Need to copy the catalogs since they are modified later.
682 inputVisitSummary = [cat.copy(deep=True) for cat in self.inputVisitSummary]
683 # Set one WCS to be None
684 testVisit = 0
685 testDetector = 20
686 inputVisitSummary[testVisit][testDetector].setWcs(None)
688 outputs = self.task.run(
689 self.inputCatalogRefs,
690 inputVisitSummary,
691 instrumentName=self.instrumentName,
692 refEpoch=self.refEpoch,
693 refObjectLoader=self.refObjectLoader,
694 )
696 # Check that the fit WCS for the extension with input WCS=None returns
697 # finite sky values.
698 testWcs = outputs.outputWcss[self.testVisits[testVisit]][testDetector].getWcs()
699 testSky = testWcs.pixelToSky(0, 0)
700 self.assertTrue(testSky.isFinite())
702 def test_inputCameraModel(self):
703 """Test running task with an input camera model, and check that true
704 object coordinates are recovered.
705 """
706 config = copy(self.config)
707 config.saveCameraModel = False
708 config.useInputCameraModel = True
709 task = GbdesAstrometricFitTask(config=config)
710 with open(os.path.join(self.datadir, "sample_camera_model.yaml"), "r") as f:
711 cameraModel = yaml.load(f, Loader=yaml.Loader)
713 outputs = task.run(
714 self.inputCatalogRefs,
715 self.inputVisitSummary,
716 instrumentName=self.instrumentName,
717 refObjectLoader=self.refObjectLoader,
718 inputCameraModel=cameraModel,
719 )
721 # Check that the output WCS is close (not necessarily exactly equal) to
722 # the result when fitting the full model.
723 for v, visit in enumerate(self.testVisits):
724 visitSummary = self.inputVisitSummary[v]
725 outputWcsCatalog = outputs.outputWcss[visit]
726 visitSources = self.inputCatalogRefs[v].get()
727 for d, detectorRow in enumerate(visitSummary):
728 detectorId = detectorRow["id"]
729 fitwcs = outputWcsCatalog[d].getWcs()
730 detSources = visitSources[visitSources["detector"] == detectorId]
731 fitRA, fitDec = fitwcs.pixelToSkyArray(detSources["x"], detSources["y"], degrees=True)
732 dRA = fitRA - detSources["trueRA"]
733 dDec = fitDec - detSources["trueDec"]
734 # Check that input coordinates match the output coordinates
735 self.assertAlmostEqual(np.mean(dRA), 0)
736 self.assertAlmostEqual(np.std(dRA), 0)
737 self.assertAlmostEqual(np.mean(dDec), 0)
738 self.assertAlmostEqual(np.std(dDec), 0)
740 def test_useColor(self):
741 """Test running task with color catalog and DCR fitting."""
742 config = copy(self.config)
743 config.useColor = True
745 task = GbdesAstrometricFitTask(config=config)
746 outputs = task.run(
747 self.inputCatalogRefs,
748 self.inputVisitSummary,
749 instrumentName=self.instrumentName,
750 refObjectLoader=self.refObjectLoader,
751 colorCatalog=self.colorCatalog,
752 )
753 self.assertEqual(len(outputs.colorParams["visit"]), len(self.testVisits))
756class TestGbdesGlobalAstrometricFit(TestGbdesAstrometricFit):
757 @classmethod
758 def setUpClass(cls):
759 # Set random seed
760 np.random.seed(1234)
762 # Fraction of simulated stars in the reference catalog and science
763 # exposures
764 inReferenceFraction = 1
765 inScienceFraction = 1
767 # Make fake data
768 cls.datadir = os.path.join(TESTDIR, "data")
770 cls.nFields = 2
771 cls.instrumentName = "HSC"
772 cls.refEpoch = 57205.5
774 # Make test inputVisitSummary. VisitSummaryTables are taken from
775 # collection HSC/runs/RC2/w_2022_20/DM-34794
776 cls.testVisits = [1176, 17900, 17930, 17934, 36434, 36460, 36494, 36446]
777 cls.inputVisitSummary = []
778 for testVisit in cls.testVisits:
779 visSum = afwTable.ExposureCatalog.readFits(
780 os.path.join(cls.datadir, f"visitSummary_{testVisit}.fits")
781 )
782 cls.inputVisitSummary.append(visSum)
784 cls.config = GbdesGlobalAstrometricFitConfig()
785 cls.config.systematicError = 0
786 cls.config.devicePolyOrder = 4
787 cls.config.exposurePolyOrder = 6
788 cls.config.fitReserveFraction = 0
789 cls.config.fitReserveRandomSeed = 1234
790 cls.config.saveModelParams = True
791 cls.config.saveCameraModel = True
792 cls.config.applyRefCatProperMotion = False
793 cls.task = GbdesGlobalAstrometricFitTask(config=cls.config)
795 cls.fields, cls.fieldRegions = cls.task._prep_sky(cls.inputVisitSummary)
797 (
798 cls.exposureInfo,
799 cls.exposuresHelper,
800 cls.extensionInfo,
801 cls.instruments,
802 ) = cls.task._get_exposure_info(cls.inputVisitSummary, fieldRegions=cls.fieldRegions)
804 refDataIds, deferredRefCats = [], []
805 allStarIds = []
806 allStarRAs = []
807 allStarDecs = []
808 for region in cls.fieldRegions.values():
809 # Bounding box of observations:
810 bbox = region.getBoundingBox()
811 raMin = bbox.getLon().getA().asDegrees()
812 raMax = bbox.getLon().getB().asDegrees()
813 decMin = bbox.getLat().getA().asDegrees()
814 decMax = bbox.getLat().getB().asDegrees()
816 # Make random set of data in a bounding box determined by input
817 # visits
818 cls.nStars, starIds, starRAs, starDecs = cls._make_simulated_stars(
819 raMin, raMax, decMin, decMax, cls.config.matchRadius
820 )
822 allStarIds.append(starIds)
823 allStarRAs.append(starRAs)
824 allStarDecs.append(starDecs)
826 corners = [
827 lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector()
828 for (ra, dec) in zip([raMin, raMax, raMax, raMin], [decMin, decMin, decMax, decMax])
829 ]
830 conv_box = lsst.sphgeom.ConvexPolygon.convexHull(corners)
831 # Make a reference catalog that will be loaded into
832 # ReferenceObjectLoader
833 refDataId, deferredRefCat = cls._make_refCat(
834 starIds, starRAs, starDecs, inReferenceFraction, conv_box
835 )
836 refDataIds.append(refDataId)
837 deferredRefCats.append(deferredRefCat)
839 cls.refObjectLoader = ReferenceObjectLoader(refDataIds, deferredRefCats)
840 cls.refObjectLoader.config.requireProperMotion = False
841 cls.refObjectLoader.config.anyFilterMapsToThis = "test_filter"
843 cls.task.refObjectLoader = cls.refObjectLoader
845 allRefObjects, allRefCovariances = {}, {}
846 for f, fieldRegion in cls.fieldRegions.items():
847 refObjects, refCovariance = cls.task._load_refcat(
848 cls.refObjectLoader, cls.extensionInfo, epoch=cls.exposureInfo.medianEpoch, region=fieldRegion
849 )
850 allRefObjects[f] = refObjects
851 allRefCovariances[f] = refCovariance
852 cls.refObjects = allRefObjects
854 # Get True WCS for stars:
855 with open(os.path.join(cls.datadir, "sample_global_wcs.yaml"), "r") as f:
856 cls.trueModel = yaml.load(f, Loader=yaml.Loader)
858 cls.trueWCSs = cls._make_wcs(cls.trueModel, cls.inputVisitSummary)
860 cls.isolatedStarCatalogs, cls.isolatedStarSources = cls._make_isolatedStars(
861 allStarIds, allStarRAs, allStarDecs, cls.trueWCSs, inScienceFraction
862 )
864 cls.colorCatalog = cls._make_colors(np.concatenate(allStarRAs), np.concatenate(allStarDecs))
866 cls.outputs = cls.task.run(
867 cls.inputVisitSummary,
868 cls.isolatedStarSources,
869 cls.isolatedStarCatalogs,
870 instrumentName=cls.instrumentName,
871 refEpoch=cls.refEpoch,
872 refObjectLoader=cls.refObjectLoader,
873 )
875 @classmethod
876 def _make_isolatedStars(cls, allStarIds, allStarRAs, allStarDecs, trueWCSs, inScienceFraction):
877 """Given a subset of the simulated data, make source catalogs and star
878 catalogs.
880 This takes the true WCSs to go from the RA and Decs of the simulated
881 stars to pixel coordinates for a given visit and detector. If those
882 pixel coordinates are within the bounding box of the detector, the
883 source and visit information is put in the corresponding catalog of
884 isolated star sources.
886 Parameters
887 ----------
888 allStarIds : `np.ndarray` [`int`]
889 Source ids for the simulated stars
890 allStarRas : `np.ndarray` [`float`]
891 RAs of the simulated stars
892 allStarDecs : `np.ndarray` [`float`]
893 Decs of the simulated stars
894 trueWCSs : `list` [`lsst.afw.geom.SkyWcs`]
895 WCS with which to simulate the source pixel coordinates
896 inReferenceFraction : float
897 Percentage of simulated stars to include in reference catalog
899 Returns
900 -------
901 isolatedStarCatalogRefs : `list`
902 [`lsst.pipe.base.InMemoryDatasetHandle`]
903 List of references to isolated star catalogs.
904 isolatedStarSourceRefs : `list`
905 [`lsst.pipe.base.InMemoryDatasetHandle`]
906 List of references to isolated star sources.
907 """
908 bbox = lsst.geom.BoxD(
909 lsst.geom.Point2D(
910 cls.inputVisitSummary[0][0]["bbox_min_x"], cls.inputVisitSummary[0][0]["bbox_min_y"]
911 ),
912 lsst.geom.Point2D(
913 cls.inputVisitSummary[0][0]["bbox_max_x"], cls.inputVisitSummary[0][0]["bbox_max_y"]
914 ),
915 )
916 bboxCorners = bbox.getCorners()
918 isolatedStarCatalogRefs = []
919 isolatedStarSourceRefs = []
920 for i in range(len(allStarIds)):
921 starIds = allStarIds[i]
922 starRAs = allStarRAs[i]
923 starDecs = allStarDecs[i]
924 isolatedStarCatalog = pd.DataFrame({"ra": starRAs, "dec": starDecs}, index=starIds)
925 isolatedStarCatalogRefs.append(
926 InMemoryDatasetHandle(isolatedStarCatalog, storageClass="DataFrame", dataId={"tract": 0})
927 )
928 sourceCats = []
929 for v, visit in enumerate(cls.testVisits):
930 nVisStars = int(cls.nStars * inScienceFraction)
931 visitStarIndices = np.random.choice(cls.nStars, nVisStars, replace=False)
932 visitStarIds = starIds[visitStarIndices]
933 visitStarRas = starRAs[visitStarIndices]
934 visitStarDecs = starDecs[visitStarIndices]
935 for detector in trueWCSs[visit]:
936 detWcs = detector.getWcs()
937 detectorId = detector["id"]
938 radecCorners = detWcs.pixelToSky(bboxCorners)
939 detectorFootprint = sphgeom.ConvexPolygon([rd.getVector() for rd in radecCorners])
940 detectorIndices = detectorFootprint.contains(
941 (visitStarRas * u.degree).to(u.radian), (visitStarDecs * u.degree).to(u.radian)
942 )
943 nDetectorStars = detectorIndices.sum()
944 detectorArray = np.ones(nDetectorStars, dtype=int) * detector["id"]
945 visitArray = np.ones(nDetectorStars, dtype=int) * visit
947 ones_like = np.ones(nDetectorStars)
948 zeros_like = np.zeros(nDetectorStars, dtype=bool)
950 x, y = detWcs.skyToPixelArray(
951 visitStarRas[detectorIndices], visitStarDecs[detectorIndices], degrees=True
952 )
954 origWcs = (cls.inputVisitSummary[v][cls.inputVisitSummary[v]["id"] == detectorId])[
955 0
956 ].getWcs()
957 inputRa, inputDec = origWcs.pixelToSkyArray(x, y, degrees=True)
959 sourceDict = {}
960 sourceDict["detector"] = detectorArray
961 sourceDict["visit"] = visitArray
962 sourceDict["obj_index"] = visitStarIds[detectorIndices]
963 sourceDict["x"] = x
964 sourceDict["y"] = y
965 sourceDict["xErr"] = 1e-3 * ones_like
966 sourceDict["yErr"] = 1e-3 * ones_like
967 sourceDict["inputRA"] = inputRa
968 sourceDict["inputDec"] = inputDec
969 sourceDict["trueRA"] = visitStarRas[detectorIndices]
970 sourceDict["trueDec"] = visitStarDecs[detectorIndices]
971 for key in ["apFlux_12_0_flux", "apFlux_12_0_instFlux", "ixx", "iyy"]:
972 sourceDict[key] = ones_like
973 for key in [
974 "pixelFlags_edge",
975 "pixelFlags_saturated",
976 "pixelFlags_interpolatedCenter",
977 "pixelFlags_interpolated",
978 "pixelFlags_crCenter",
979 "pixelFlags_bad",
980 "hsmPsfMoments_flag",
981 "apFlux_12_0_flag",
982 "extendedness",
983 "parentSourceId",
984 "deblend_nChild",
985 "ixy",
986 ]:
987 sourceDict[key] = zeros_like
988 sourceDict["apFlux_12_0_instFluxErr"] = 1e-2 * ones_like
989 sourceDict["detect_isPrimary"] = ones_like.astype(bool)
991 sourceCat = pd.DataFrame(sourceDict)
992 sourceCats.append(sourceCat)
994 isolatedStarSourceTable = pd.concat(sourceCats, ignore_index=True)
995 isolatedStarSourceTable = isolatedStarSourceTable.sort_values(by=["obj_index"])
996 isolatedStarSourceRefs.append(
997 InMemoryDatasetHandle(isolatedStarSourceTable, storageClass="DataFrame", dataId={"tract": 0})
998 )
1000 return isolatedStarCatalogRefs, isolatedStarSourceRefs
1002 def test_loading_and_association(self):
1003 """Test that associated objects actually correspond to the same
1004 simulated object."""
1005 associations, sourceDict = self.task._associate_from_isolated_sources(
1006 self.isolatedStarSources, self.isolatedStarCatalogs, self.extensionInfo, self.refObjects
1007 )
1009 object_groups = np.flatnonzero(np.array(associations.sequence) == 0)
1010 for i in range(len(object_groups) - 1)[:10]:
1011 ras, decs = [], []
1012 for ind in np.arange(object_groups[i], object_groups[i + 1]):
1013 visit = self.extensionInfo.visit[associations.extn[ind]]
1014 detectorInd = self.extensionInfo.detectorIndex[associations.extn[ind]]
1015 detector = self.extensionInfo.detector[associations.extn[ind]]
1016 if detectorInd == -1:
1017 ra = self.refObjects[visit * -1]["ra"][associations.obj[ind]]
1018 dec = self.refObjects[visit * -1]["dec"][associations.obj[ind]]
1019 ras.append(ra)
1020 decs.append(dec)
1021 else:
1022 x = sourceDict[visit][detector]["x"][associations.obj[ind]]
1023 y = sourceDict[visit][detector]["y"][associations.obj[ind]]
1024 ra, dec = self.trueWCSs[visit][detectorInd].getWcs().pixelToSky(x, y)
1025 ras.append(ra.asDegrees())
1026 decs.append(dec.asDegrees())
1027 np.testing.assert_allclose(ras, ras[0])
1028 np.testing.assert_allclose(decs, decs[0])
1030 def test_refCatLoader(self):
1031 """Test loading objects from the refCat in each of the fields."""
1033 for region in self.fieldRegions.values():
1034 refCat, refCov = self.task._load_refcat(
1035 self.refObjectLoader, self.extensionInfo, region=region, epoch=self.exposureInfo.medianEpoch
1036 )
1037 assert len(refCat) > 0
1039 def test_make_outputs(self):
1040 """Test that the run method recovers the input model parameters."""
1041 for isolatedStarSourceRef in self.isolatedStarSources:
1042 iss = isolatedStarSourceRef.get()
1043 visits = np.unique(iss["visit"])
1044 for v, visit in enumerate(visits):
1045 outputWcsCatalog = self.outputs.outputWcss[visit]
1046 visitSources = iss[iss["visit"] == visit]
1047 detectors = outputWcsCatalog["id"]
1048 for d, detectorId in enumerate(detectors):
1049 fitwcs = outputWcsCatalog[d].getWcs()
1050 detSources = visitSources[visitSources["detector"] == detectorId]
1051 fitRA, fitDec = fitwcs.pixelToSkyArray(detSources["x"], detSources["y"], degrees=True)
1052 dRA = fitRA - detSources["trueRA"]
1053 dDec = fitDec - detSources["trueDec"]
1054 # Check that input coordinates match the output coordinates
1055 self.assertAlmostEqual(np.mean(dRA), 0)
1056 self.assertAlmostEqual(np.std(dRA), 0)
1057 self.assertAlmostEqual(np.mean(dDec), 0)
1058 self.assertAlmostEqual(np.std(dDec), 0)
1060 def test_missingWcs(self):
1061 """Test that task does not fail when the input WCS is None for one
1062 extension and that the fit WCS for that extension returns a finite
1063 result.
1064 """
1065 # Need to copy the catalogs since they are modified later.
1066 inputVisitSummary = [cat.copy(deep=True) for cat in self.inputVisitSummary]
1067 # Set one WCS to be None
1068 testVisit = 0
1069 testDetector = 20
1070 inputVisitSummary[testVisit][testDetector].setWcs(None)
1072 outputs = self.task.run(
1073 inputVisitSummary,
1074 self.isolatedStarSources,
1075 self.isolatedStarCatalogs,
1076 instrumentName=self.instrumentName,
1077 refEpoch=self.refEpoch,
1078 refObjectLoader=self.refObjectLoader,
1079 )
1081 # Check that the fit WCS for the extension with input WCS=None returns
1082 # finite sky values.
1083 testWcs = outputs.outputWcss[self.testVisits[testVisit]][testDetector].getWcs()
1084 testSky = testWcs.pixelToSky(0, 0)
1085 self.assertTrue(testSky.isFinite())
1087 def test_inputCameraModel(self):
1088 """Test running task with an input camera model, and check that true
1089 object coordinates are recovered.
1090 """
1091 config = copy(self.config)
1092 config.saveCameraModel = False
1093 config.useInputCameraModel = True
1094 task = GbdesGlobalAstrometricFitTask(config=config)
1095 with open(os.path.join(self.datadir, "sample_global_camera_model.yaml"), "r") as f:
1096 cameraModel = yaml.load(f, Loader=yaml.Loader)
1098 outputs = task.run(
1099 self.inputVisitSummary,
1100 self.isolatedStarSources,
1101 self.isolatedStarCatalogs,
1102 instrumentName=self.instrumentName,
1103 refObjectLoader=self.refObjectLoader,
1104 inputCameraModel=cameraModel,
1105 )
1107 # Check that the output WCS is close (not necessarily exactly equal) to
1108 # the result when fitting the full model.
1109 for isolatedStarSourceRef in self.isolatedStarSources:
1110 iss = isolatedStarSourceRef.get()
1111 visits = np.unique(iss["visit"])
1112 for v, visit in enumerate(visits):
1113 outputWcsCatalog = outputs.outputWcss[visit]
1114 visitSources = iss[iss["visit"] == visit]
1115 detectors = outputWcsCatalog["id"]
1116 for d, detectorId in enumerate(detectors):
1117 fitwcs = outputWcsCatalog[d].getWcs()
1118 detSources = visitSources[visitSources["detector"] == detectorId]
1119 fitRA, fitDec = fitwcs.pixelToSkyArray(detSources["x"], detSources["y"], degrees=True)
1120 dRA = fitRA - detSources["trueRA"]
1121 dDec = fitDec - detSources["trueDec"]
1122 # Check that input coordinates match the output coordinates
1123 self.assertAlmostEqual(np.mean(dRA), 0, places=6)
1124 self.assertAlmostEqual(np.std(dRA), 0)
1125 self.assertAlmostEqual(np.mean(dDec), 0, places=6)
1126 self.assertAlmostEqual(np.std(dDec), 0)
1128 def test_useColor(self):
1129 """Test running task with color catalog and DCR fitting."""
1130 config = copy(self.config)
1131 config.useColor = True
1133 task = GbdesGlobalAstrometricFitTask(config=config)
1135 outputs = task.run(
1136 self.inputVisitSummary,
1137 self.isolatedStarSources,
1138 self.isolatedStarCatalogs,
1139 instrumentName=self.instrumentName,
1140 refObjectLoader=self.refObjectLoader,
1141 colorCatalog=self.colorCatalog,
1142 )
1143 self.assertEqual(len(outputs.colorParams["visit"]), len(self.testVisits))
1146class TestApparentMotion(lsst.utils.tests.TestCase):
1147 """This class simulates an array of objects with random proper motions and
1148 parallaxes and compares the results of
1149 `lsst.drp.tasks.gbdesAstrometricFit.calculate_apparent_motion` against the
1150 results calculated by astropy."""
1152 def setUp(self):
1153 np.random.seed(12345)
1155 self.ras = (np.random.rand(10) + 150.0) * u.degree
1156 self.decs = (np.random.rand(10) + 2.5) * u.degree
1158 self.pmRAs = (np.random.random(10) * 10) * u.mas / u.yr
1159 self.pmDecs = (np.random.random(10) * 10) * u.mas / u.yr
1161 self.parallaxes = (abs(np.random.random(10) * 5)) * u.mas
1163 self.mjds = astropy.time.Time(np.random.rand(10) * 20 + 57000, format="mjd", scale="tai")
1165 self.refEpoch = astropy.time.Time(57100, format="mjd", scale="tai")
1167 def test_proper_motion(self):
1168 """Calculate the change in position due to proper motion only for a
1169 given time change, and compare the results against astropy.
1170 """
1171 # Set parallaxes to zero to test proper motion part by itself.
1172 parallaxes = np.zeros(len(self.ras)) * u.mas
1174 data = astropy.table.Table(
1175 {
1176 "ra": self.ras,
1177 "dec": self.decs,
1178 "pmRA": self.pmRAs,
1179 "pmDec": self.pmDecs,
1180 "parallax": parallaxes,
1181 "MJD": self.mjds,
1182 }
1183 )
1184 raCorrection, decCorrection = calculate_apparent_motion(data, self.refEpoch)
1186 pmRACosDec = self.pmRAs * np.cos(self.decs)
1187 coords = SkyCoord(
1188 ra=self.ras,
1189 dec=self.decs,
1190 pm_ra_cosdec=pmRACosDec,
1191 pm_dec=self.pmDecs,
1192 obstime=self.refEpoch,
1193 )
1195 newCoords = coords.apply_space_motion(new_obstime=self.mjds)
1196 astropyDRA = newCoords.ra.degree - coords.ra.degree
1197 astropyDDec = newCoords.dec.degree - coords.dec.degree
1199 # The proper motion values are on the order of 0-3 mas/yr and the
1200 # differences are on the scale of 1e-8 mas.
1201 self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=1e-6)
1202 self.assertFloatsAlmostEqual(astropyDDec, decCorrection.value, rtol=1e-6)
1204 def test_parallax(self):
1205 """Calculate the change in position due to parallax for a given time
1206 change and compare the results against astropy. Astropy will not give
1207 the parallax part separate from other sources of space motion, such as
1208 annual aberration, but it can be obtained by comparing the calculated
1209 positions of an object with a given parallax and one with a much
1210 further distance. The result is still only approximately the same, but
1211 is close enough for accuracy needed here.
1212 """
1213 # Set proper motions to zero to test parallax correction by itself.
1214 pms = np.zeros(len(self.ras)) * u.mas / u.yr
1216 data = astropy.table.Table(
1217 {
1218 "ra": self.ras,
1219 "dec": self.decs,
1220 "pmRA": pms,
1221 "pmDec": pms,
1222 "parallax": self.parallaxes,
1223 "MJD": self.mjds,
1224 }
1225 )
1226 raCorrection, decCorrection = calculate_apparent_motion(data, self.refEpoch)
1228 coords = SkyCoord(
1229 ra=self.ras,
1230 dec=self.decs,
1231 pm_ra_cosdec=pms,
1232 pm_dec=pms,
1233 distance=Distance(parallax=self.parallaxes),
1234 obstime=self.refEpoch,
1235 )
1236 # Make another set of objects at the same coordinates but a much
1237 # greater distant in order to isolate parallax when applying space
1238 # motion.
1239 refCoords = SkyCoord(
1240 ra=self.ras,
1241 dec=self.decs,
1242 pm_ra_cosdec=pms,
1243 pm_dec=pms,
1244 distance=1e10 * Distance(parallax=self.parallaxes),
1245 obstime=self.refEpoch,
1246 )
1248 newCoords = coords.apply_space_motion(new_obstime=self.mjds).transform_to("gcrs")
1249 refNewCoords = refCoords.apply_space_motion(new_obstime=self.mjds).transform_to("gcrs")
1250 astropyDRA = newCoords.ra.degree - refNewCoords.ra.degree
1251 astropyDDec = newCoords.dec.degree - refNewCoords.dec.degree
1253 # The parallax values are on the scale of 1-3 mas, and the difference
1254 # between the two calculation methods is ~0.05 mas.
1255 self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=2e-2)
1256 self.assertFloatsAlmostEqual(astropyDDec, decCorrection.value, rtol=2e-2)
1259def setup_module(module):
1260 lsst.utils.tests.init()
1263if __name__ == "__main__": 1263 ↛ 1264line 1263 didn't jump to line 1264 because the condition on line 1263 was never true
1264 lsst.utils.tests.init()
1265 unittest.main()