Coverage for tests/test_gbdesAstrometricFit.py: 9%
602 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-18 00:46 -0800
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.task = GbdesAstrometricFitTask(config=cls.config)
95 (
96 cls.exposureInfo,
97 cls.exposuresHelper,
98 cls.extensionInfo,
99 cls.instruments,
100 ) = cls.task._get_exposure_info(cls.inputVisitSummary, refEpoch=cls.refEpoch)
102 cls.fields, cls.fieldCenter, cls.fieldRadius = cls.task._prep_sky(
103 cls.inputVisitSummary, cls.exposureInfo.medianEpoch
104 )
106 # Bounding box of observations:
107 raMins, raMaxs = [], []
108 decMins, decMaxs = [], []
109 for visSum in cls.inputVisitSummary:
110 raMins.append(visSum["raCorners"].min())
111 raMaxs.append(visSum["raCorners"].max())
112 decMins.append(visSum["decCorners"].min())
113 decMaxs.append(visSum["decCorners"].max())
114 raMin = min(raMins)
115 raMax = max(raMaxs)
116 decMin = min(decMins)
117 decMax = max(decMaxs)
119 corners = [
120 lsst.geom.SpherePoint(raMin, decMin, lsst.geom.degrees).getVector(),
121 lsst.geom.SpherePoint(raMax, decMin, lsst.geom.degrees).getVector(),
122 lsst.geom.SpherePoint(raMax, decMax, lsst.geom.degrees).getVector(),
123 lsst.geom.SpherePoint(raMin, decMax, lsst.geom.degrees).getVector(),
124 ]
125 cls.boundingPolygon = sphgeom.ConvexPolygon(corners)
127 # Make random set of data in a bounding box determined by input visits
128 # Make wcs objects for the "true" model
129 cls.nStars, starIds, starRAs, starDecs = cls._make_simulated_stars(
130 raMin, raMax, decMin, decMax, cls.config.matchRadius
131 )
133 # Fraction of simulated stars in the reference catalog and science
134 # exposures
135 inReferenceFraction = 1
136 inScienceFraction = 1
138 # Make a reference catalog and load it into ReferenceObjectLoader
139 refDataId, deferredRefCat = cls._make_refCat(
140 starIds, starRAs, starDecs, inReferenceFraction, cls.boundingPolygon
141 )
142 cls.refObjectLoader = ReferenceObjectLoader([refDataId], [deferredRefCat])
143 cls.refObjectLoader.config.requireProperMotion = False
144 cls.refObjectLoader.config.anyFilterMapsToThis = "test_filter"
146 cls.task.refObjectLoader = cls.refObjectLoader
148 # Get True WCS for stars:
149 with open(os.path.join(cls.datadir, "sample_wcs.yaml"), "r") as f:
150 cls.trueModel = yaml.load(f, Loader=yaml.Loader)
152 trueWCSs = cls._make_wcs(cls.trueModel, cls.inputVisitSummary)
154 # Make source catalogs:
155 cls.inputCatalogRefs = cls._make_sourceCat(starIds, starRAs, starDecs, trueWCSs, inScienceFraction)
157 cls.colorCatalog = cls._make_colors(starRAs, starDecs)
159 cls.outputs = cls.task.run(
160 cls.inputCatalogRefs,
161 cls.inputVisitSummary,
162 instrumentName=cls.instrumentName,
163 refEpoch=cls.refEpoch,
164 refObjectLoader=cls.refObjectLoader,
165 )
167 @staticmethod
168 def _make_simulated_stars(raMin, raMax, decMin, decMax, matchRadius, nStars=10000):
169 """Generate random positions for "stars" in an RA/Dec box.
171 Parameters
172 ----------
173 raMin : `float`
174 Minimum RA for simulated stars.
175 raMax : `float`
176 Maximum RA for simulated stars.
177 decMin : `float`
178 Minimum Dec for simulated stars.
179 decMax : `float`
180 Maximum Dec for simulated stars.
181 matchRadius : `float`
182 Minimum allowed distance in arcsec between stars.
183 nStars : `int`, optional
184 Number of stars to simulate. Final number will be lower if any
185 too-close stars are dropped.
187 Returns
188 -------
189 nStars : `int`
190 Number of stars simulated.
191 starIds: `np.ndarray`
192 Unique identification number for stars.
193 starRAs: `np.ndarray`
194 Simulated Right Ascensions.
195 starDecs: `np.ndarray`
196 Simulated Declination.
197 """
198 starIds = np.arange(nStars)
199 starRAs = np.random.random(nStars) * (raMax - raMin) + raMin
200 starDecs = np.random.random(nStars) * (decMax - decMin) + decMin
201 # Remove neighbors:
202 with Matcher(starRAs, starDecs) as matcher:
203 idx = matcher.query_groups(matchRadius / 3600.0, min_match=2)
204 if len(idx) > 0:
205 neighbors = np.unique(np.concatenate(idx))
206 starRAs = np.delete(starRAs, neighbors)
207 starDecs = np.delete(starDecs, neighbors)
208 nStars = len(starRAs)
209 starIds = np.arange(nStars)
211 return nStars, starIds, starRAs, starDecs
213 @classmethod
214 def _make_refCat(cls, starIds, starRas, starDecs, inReferenceFraction, bounds):
215 """Make reference catalog from a subset of the simulated data
217 Parameters
218 ----------
219 starIds : `np.ndarray` [`int`]
220 Source ids for the simulated stars
221 starRas : `np.ndarray` [`float`]
222 RAs of the simulated stars
223 starDecs : `np.ndarray` [`float`]
224 Decs of the simulated stars
225 inReferenceFraction : float
226 Percentage of simulated stars to include in reference catalog
227 bounds : `lsst.sphgeom.ConvexPolygon`
228 Boundary of the reference catalog region
230 Returns
231 -------
232 refDataId : `lsst.meas.algorithms.testUtils.MockRefcatDataId`
233 Object that replicates the functionality of a dataId.
234 deferredRefCat : `lsst.pipe.base.InMemoryDatasetHandle`
235 Dataset handle for reference catalog.
236 """
237 nRefs = int(cls.nStars * inReferenceFraction)
238 refStarIndices = np.random.choice(cls.nStars, nRefs, replace=False)
239 # Make simpleCatalog to hold data, create datasetRef with `region`
240 # determined by bounding box used in above simulate.
241 refSchema = afwTable.SimpleTable.makeMinimalSchema()
242 idKey = refSchema.addField("sourceId", type="I")
243 fluxKey = refSchema.addField("test_filter_flux", units="nJy", type=np.float64)
244 raErrKey = refSchema.addField("coord_raErr", units="rad", type=np.float64)
245 decErrKey = refSchema.addField("coord_decErr", units="rad", type=np.float64)
246 pmraErrKey = refSchema.addField("pm_raErr", units="rad2 / yr", type=np.float64)
247 pmdecErrKey = refSchema.addField("pm_decErr", units="rad2 / yr", type=np.float64)
248 refCat = afwTable.SimpleCatalog(refSchema)
249 ref_md = PropertyList()
250 ref_md.set("REFCAT_FORMAT_VERSION", 1)
251 refCat.table.setMetadata(ref_md)
252 for i in refStarIndices:
253 record = refCat.addNew()
254 record.set(idKey, starIds[i])
255 record.setRa(lsst.geom.Angle(starRas[i], lsst.geom.degrees))
256 record.setDec(lsst.geom.Angle(starDecs[i], lsst.geom.degrees))
257 record.set(fluxKey, 1)
258 record.set(raErrKey, 0.00001)
259 record.set(decErrKey, 0.00001)
260 record.set(pmraErrKey, 1e-9)
261 record.set(pmdecErrKey, 1e-9)
262 refDataId = MockRefcatDataId(bounds)
263 deferredRefCat = InMemoryDatasetHandle(refCat, storageClass="SourceCatalog", htm7="mockRefCat")
265 return refDataId, deferredRefCat
267 @classmethod
268 def _make_sourceCat(cls, starIds, starRas, starDecs, trueWCSs, inScienceFraction):
269 """Make a `pd.DataFrame` catalog with the columns needed for the
270 object selector.
272 Parameters
273 ----------
274 starIds : `np.ndarray` [`int`]
275 Source ids for the simulated stars
276 starRas : `np.ndarray` [`float`]
277 RAs of the simulated stars
278 starDecs : `np.ndarray` [`float`]
279 Decs of the simulated stars
280 trueWCSs : `list` [`lsst.afw.geom.SkyWcs`]
281 WCS with which to simulate the source pixel coordinates
282 inReferenceFraction : float
283 Percentage of simulated stars to include in reference catalog
285 Returns
286 -------
287 sourceCat : `list` [`lsst.pipe.base.InMemoryDatasetHandle`]
288 List of reference to source catalogs.
289 """
290 inputCatalogRefs = []
291 # Take a subset of the simulated data
292 # Use true wcs objects to put simulated data into ccds
293 bbox = lsst.geom.BoxD(
294 lsst.geom.Point2D(
295 cls.inputVisitSummary[0][0]["bbox_min_x"], cls.inputVisitSummary[0][0]["bbox_min_y"]
296 ),
297 lsst.geom.Point2D(
298 cls.inputVisitSummary[0][0]["bbox_max_x"], cls.inputVisitSummary[0][0]["bbox_max_y"]
299 ),
300 )
301 bboxCorners = bbox.getCorners()
302 cls.inputCatalogRefs = []
303 for v, visit in enumerate(cls.testVisits):
304 nVisStars = int(cls.nStars * inScienceFraction)
305 visitStarIndices = np.random.choice(cls.nStars, nVisStars, replace=False)
306 visitStarIds = starIds[visitStarIndices]
307 visitStarRas = starRas[visitStarIndices]
308 visitStarDecs = starDecs[visitStarIndices]
309 sourceCats = []
310 for detector in trueWCSs[visit]:
311 detWcs = detector.getWcs()
312 detectorId = detector["id"]
313 radecCorners = detWcs.pixelToSky(bboxCorners)
314 detectorFootprint = sphgeom.ConvexPolygon([rd.getVector() for rd in radecCorners])
315 detectorIndices = detectorFootprint.contains(
316 (visitStarRas * u.degree).to(u.radian), (visitStarDecs * u.degree).to(u.radian)
317 )
318 nDetectorStars = detectorIndices.sum()
319 detectorArray = np.ones(nDetectorStars, dtype=bool) * detector["id"]
321 ones_like = np.ones(nDetectorStars)
322 zeros_like = np.zeros(nDetectorStars, dtype=bool)
324 x, y = detWcs.skyToPixelArray(
325 visitStarRas[detectorIndices], visitStarDecs[detectorIndices], degrees=True
326 )
328 origWcs = (cls.inputVisitSummary[v][cls.inputVisitSummary[v]["id"] == detectorId])[0].getWcs()
329 inputRa, inputDec = origWcs.pixelToSkyArray(x, y, degrees=True)
331 sourceDict = {}
332 sourceDict["detector"] = detectorArray
333 sourceDict["sourceId"] = visitStarIds[detectorIndices]
334 sourceDict["x"] = x
335 sourceDict["y"] = y
336 sourceDict["xErr"] = 1e-3 * ones_like
337 sourceDict["yErr"] = 1e-3 * ones_like
338 sourceDict["inputRA"] = inputRa
339 sourceDict["inputDec"] = inputDec
340 sourceDict["trueRA"] = visitStarRas[detectorIndices]
341 sourceDict["trueDec"] = visitStarDecs[detectorIndices]
342 for key in ["apFlux_12_0_flux", "apFlux_12_0_instFlux", "ixx", "iyy"]:
343 sourceDict[key] = ones_like
344 for key in [
345 "pixelFlags_edge",
346 "pixelFlags_saturated",
347 "pixelFlags_interpolatedCenter",
348 "pixelFlags_interpolated",
349 "pixelFlags_crCenter",
350 "pixelFlags_bad",
351 "hsmPsfMoments_flag",
352 "apFlux_12_0_flag",
353 "extendedness",
354 "sizeExtendedness",
355 "parentSourceId",
356 "deblend_nChild",
357 "ixy",
358 ]:
359 sourceDict[key] = zeros_like
360 sourceDict["apFlux_12_0_instFluxErr"] = 1e-2 * ones_like
361 sourceDict["detect_isPrimary"] = ones_like.astype(bool)
363 sourceCat = pd.DataFrame(sourceDict)
364 sourceCats.append(sourceCat)
366 visitSourceTable = pd.concat(sourceCats)
368 inputCatalogRef = InMemoryDatasetHandle(
369 visitSourceTable, storageClass="DataFrame", dataId={"visit": visit}
370 )
372 inputCatalogRefs.append(inputCatalogRef)
374 return inputCatalogRefs
376 @classmethod
377 def _make_wcs(cls, model, inputVisitSummaries):
378 """Make a `lsst.afw.geom.SkyWcs` from given model parameters
380 Parameters
381 ----------
382 model : `dict`
383 Dictionary with WCS model parameters
384 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
385 Visit summary catalogs
386 Returns
387 -------
388 catalogs : `dict` [`int`, `lsst.afw.table.ExposureCatalog`]
389 Visit summary catalogs with WCS set to input model
390 """
392 # Pixels will need to be rescaled before going into the mappings
393 xscale = inputVisitSummaries[0][0]["bbox_max_x"] - inputVisitSummaries[0][0]["bbox_min_x"]
394 yscale = inputVisitSummaries[0][0]["bbox_max_y"] - inputVisitSummaries[0][0]["bbox_min_y"]
396 catalogs = {}
397 schema = lsst.afw.table.ExposureTable.makeMinimalSchema()
398 schema.addField("visit", type="L", doc="Visit number")
399 for visitSum in inputVisitSummaries:
400 visit = visitSum[0]["visit"]
401 visitMapName = f"{visit}/poly"
402 visitModel = model[visitMapName]
404 catalog = lsst.afw.table.ExposureCatalog(schema)
405 catalog.resize(len(visitSum))
406 catalog["visit"] = visit
408 raDec = visitSum[0].getVisitInfo().getBoresightRaDec()
410 visitMapType = visitModel["Type"]
411 visitDict = {"Type": visitMapType}
412 if visitMapType == "Poly":
413 mapCoefficients = visitModel["XPoly"]["Coefficients"] + visitModel["YPoly"]["Coefficients"]
414 visitDict["Coefficients"] = mapCoefficients
416 for d, detector in enumerate(visitSum):
417 detectorId = detector["id"]
418 filterName = detector["physical_filter"]
419 detectorMapName = f"{filterName}/{detectorId}/poly"
420 detectorModel = model[detectorMapName]
422 detectorMapType = detectorModel["Type"]
423 mapDict = {detectorMapName: {"Type": detectorMapType}, visitMapName: visitDict}
424 if detectorMapType == "Poly":
425 mapCoefficients = (
426 detectorModel["XPoly"]["Coefficients"] + detectorModel["YPoly"]["Coefficients"]
427 )
428 mapDict[detectorMapName]["Coefficients"] = mapCoefficients
430 outWCS = cls.task._make_afw_wcs(
431 mapDict,
432 raDec.getRa(),
433 raDec.getDec(),
434 doNormalizePixels=True,
435 xScale=xscale,
436 yScale=yscale,
437 )
438 catalog[d].setId(detectorId)
439 catalog[d].setWcs(outWCS)
441 catalog.sort()
442 catalogs[visit] = catalog
444 return catalogs
446 @classmethod
447 def _make_colors(cls, starRas, starDecs):
448 """Make a catalog with the star magnitudes.
450 Parameters
451 ----------
452 starRas : `np.ndarray` [`float`]
453 RAs of the simulated stars
454 starDecs : `np.ndarray` [`float`]
455 Decs of the simulated stars
457 Returns
458 -------
459 colorCatalog : `lsst.afw.table.SimpleCatalog`
460 Catalog with star magnitudes.
461 """
462 bands = ["g", "r", "i", "z", "y"]
463 nStars = len(starRas)
465 # Make a catalog following what is done in `fgcmCal`.
466 schema = afwTable.SimpleTable.makeMinimalSchema()
467 schema.addField(
468 "mag_std_noabs",
469 type="ArrayF",
470 doc="Standard magnitude (no absolute calibration)",
471 size=len(bands),
472 )
473 colorCatalog = afwTable.SimpleCatalog(schema)
474 colorCatalog.resize(len(starRas))
475 colorCatalog["coord_ra"] = (starRas * u.degree + 10 * np.random.randn(nStars) * u.mas).to(u.radian)
476 colorCatalog["coord_dec"] = (starDecs * u.degree + 10 * np.random.randn(nStars) * u.mas).to(u.radian)
478 magMin = 19
479 magMax = 23
480 for i in range(len(bands)):
481 colorCatalog["mag_std_noabs"][:, i] = np.random.random(nStars) * (magMax - magMin) + magMin
483 md = PropertyList()
484 md.set("BANDS", bands)
485 colorCatalog.setMetadata(md)
486 return colorCatalog
488 def test_get_exposure_info(self):
489 """Test that information for input exposures is as expected and that
490 the WCS in the class object gives approximately the same results as the
491 input `lsst.afw.geom.SkyWcs`.
492 """
494 # The total number of extensions is the number of detectors for each
495 # visit plus one for the reference catalog for each field.
496 totalExtensions = sum([len(visSum) for visSum in self.inputVisitSummary]) + self.nFields
498 self.assertEqual(totalExtensions, len(self.extensionInfo.visit))
500 taskVisits = set(self.extensionInfo.visit)
501 refVisits = np.arange(0, -1 * self.nFields, -1).tolist()
502 self.assertEqual(taskVisits, set(self.testVisits + refVisits))
504 xx = np.linspace(0, 2000, 3)
505 yy = np.linspace(0, 4000, 6)
506 xgrid, ygrid = np.meshgrid(xx, yy)
507 for visSum in self.inputVisitSummary:
508 visit = visSum[0]["visit"]
509 for detectorInfo in visSum:
510 detector = detectorInfo["id"]
511 extensionIndex = np.flatnonzero(
512 (self.extensionInfo.visit == visit) & (self.extensionInfo.detector == detector)
513 )[0]
514 fitWcs = self.extensionInfo.wcs[extensionIndex]
515 calexpWcs = detectorInfo.getWcs()
517 tanPlaneXY = np.array([fitWcs.toWorld(x, y) for (x, y) in zip(xgrid.ravel(), ygrid.ravel())])
519 calexpra, calexpdec = calexpWcs.pixelToSkyArray(xgrid.ravel(), ygrid.ravel(), degrees=True)
521 # The pixel origin maps to a position slightly off from the
522 # tangent plane origin, so we want to use this value for the
523 # tangent-plane-to-sky part of the mapping.
524 tangentPlaneCenter = fitWcs.toWorld(
525 calexpWcs.getPixelOrigin().getX(), calexpWcs.getPixelOrigin().getY()
526 )
527 tangentPlaneOrigin = lsst.geom.Point2D(tangentPlaneCenter)
528 skyOrigin = calexpWcs.pixelToSky(
529 calexpWcs.getPixelOrigin().getX(), calexpWcs.getPixelOrigin().getY()
530 )
531 cdMatrix = afwgeom.makeCdMatrix(1.0 * lsst.geom.degrees, 0.0 * lsst.geom.degrees, True)
532 iwcToSkyWcs = afwgeom.makeSkyWcs(tangentPlaneOrigin, skyOrigin, cdMatrix)
533 newRAdeg, newDecdeg = iwcToSkyWcs.pixelToSkyArray(
534 tanPlaneXY[:, 0], tanPlaneXY[:, 1], degrees=True
535 )
537 np.testing.assert_allclose(calexpra, newRAdeg)
538 np.testing.assert_allclose(calexpdec, newDecdeg)
540 def test_refCatLoader(self):
541 """Test that we can load objects from refCat"""
543 tmpAssociations = wcsfit.FoFClass(
544 self.fields,
545 self.instruments,
546 self.exposuresHelper,
547 [self.fieldRadius.asDegrees()],
548 (self.task.config.matchRadius * u.arcsec).to(u.degree).value,
549 )
551 self.task._load_refcat(
552 self.refObjectLoader,
553 self.extensionInfo,
554 associations=tmpAssociations,
555 center=self.fieldCenter,
556 radius=self.fieldRadius,
557 epoch=self.exposureInfo.medianEpoch,
558 )
560 # We have only loaded one catalog, so getting the 'matches' should just
561 # return the same objects we put in, except some random objects that
562 # are too close together.
563 tmpAssociations.sortMatches(self.fieldNumber, minMatches=1)
565 nMatches = (np.array(tmpAssociations.sequence) == 0).sum()
567 self.assertLessEqual(nMatches, self.nStars)
568 self.assertGreater(nMatches, self.nStars * 0.9)
570 def test_loading_and_association(self):
571 """Test that objects can be loaded and correctly associated."""
572 # Running `_load_catalogs_and_associate` changes the input WCSs, so
573 # recalculate them here so that the variables shared among tests are
574 # not affected.
575 _, exposuresHelper, extensionInfo, instruments = self.task._get_exposure_info(
576 self.inputVisitSummary, refEpoch=self.refEpoch
577 )
579 tmpAssociations = wcsfit.FoFClass(
580 self.fields,
581 self.instruments,
582 exposuresHelper,
583 [self.fieldRadius.asDegrees()],
584 (self.task.config.matchRadius * u.arcsec).to(u.degree).value,
585 )
586 self.task._load_catalogs_and_associate(tmpAssociations, self.inputCatalogRefs, extensionInfo)
588 tmpAssociations.sortMatches(self.fieldNumber, minMatches=2)
590 matchIds = []
591 correctMatches = []
592 for s, e, o in zip(tmpAssociations.sequence, tmpAssociations.extn, tmpAssociations.obj):
593 objVisitInd = extensionInfo.visitIndex[e]
594 objDet = extensionInfo.detector[e]
595 extnInds = self.inputCatalogRefs[objVisitInd].get()["detector"] == objDet
596 objInfo = self.inputCatalogRefs[objVisitInd].get()[extnInds].iloc[o]
597 if s == 0:
598 if len(matchIds) > 0:
599 correctMatches.append(len(set(matchIds)) == 1)
600 matchIds = []
602 matchIds.append(objInfo["sourceId"])
604 # A few matches may incorrectly associate sources because of the random
605 # positions
606 self.assertGreater(sum(correctMatches), len(correctMatches) * 0.95)
608 def test_make_outputs(self):
609 """Test that the run method recovers the input model parameters."""
610 for v, visit in enumerate(self.testVisits):
611 visitSummary = self.inputVisitSummary[v]
612 outputWcsCatalog = self.outputs.outputWcss[visit]
613 visitSources = self.inputCatalogRefs[v].get()
614 for d, detectorRow in enumerate(visitSummary):
615 detectorId = detectorRow["id"]
616 fitwcs = outputWcsCatalog[d].getWcs()
617 detSources = visitSources[visitSources["detector"] == detectorId]
618 fitRA, fitDec = fitwcs.pixelToSkyArray(detSources["x"], detSources["y"], degrees=True)
619 dRA = fitRA - detSources["trueRA"]
620 dDec = fitDec - detSources["trueDec"]
621 # Check that input coordinates match the output coordinates
622 self.assertAlmostEqual(np.mean(dRA), 0)
623 self.assertAlmostEqual(np.std(dRA), 0)
624 self.assertAlmostEqual(np.mean(dDec), 0)
625 self.assertAlmostEqual(np.std(dDec), 0)
627 def test_compute_model_params(self):
628 """Test the optional model parameters and covariance output."""
629 modelParams = pd.DataFrame(self.outputs.modelParams)
630 # Check that DataFrame is the expected size.
631 shape = modelParams.shape
632 self.assertEqual(shape[0] + 4, shape[1])
633 # Check that covariance matrix is symmetric.
634 covariance = (modelParams.iloc[:, 4:]).to_numpy()
635 np.testing.assert_allclose(covariance, covariance.T, atol=1e-18)
637 def test_run(self):
638 """Test that run method recovers the input model parameters"""
639 outputMaps = self.outputs.fitModel.mapCollection.getParamDict()
641 for v, visit in enumerate(self.testVisits):
642 visitSummary = self.inputVisitSummary[v]
643 visitMapName = f"{visit}/poly"
645 origModel = self.trueModel[visitMapName]
646 if origModel["Type"] != "Identity":
647 fitModel = outputMaps[visitMapName]
648 origXPoly = origModel["XPoly"]["Coefficients"]
649 origYPoly = origModel["YPoly"]["Coefficients"]
650 fitXPoly = fitModel[: len(origXPoly)]
651 fitYPoly = fitModel[len(origXPoly) :]
653 absDiffX = abs(fitXPoly - origXPoly)
654 absDiffY = abs(fitYPoly - origYPoly)
655 # Check that input visit model matches fit
656 np.testing.assert_array_less(absDiffX, 1e-6)
657 np.testing.assert_array_less(absDiffY, 1e-6)
658 for d, detectorRow in enumerate(visitSummary):
659 detectorId = detectorRow["id"]
660 filterName = detectorRow["physical_filter"]
661 detectorMapName = f"{filterName}/{detectorId}/poly"
662 origModel = self.trueModel[detectorMapName]
663 if (origModel["Type"] != "Identity") and (v == 0):
664 fitModel = outputMaps[detectorMapName]
665 origXPoly = origModel["XPoly"]["Coefficients"]
666 origYPoly = origModel["YPoly"]["Coefficients"]
667 fitXPoly = fitModel[: len(origXPoly)]
668 fitYPoly = fitModel[len(origXPoly) :]
669 absDiffX = abs(fitXPoly - origXPoly)
670 absDiffY = abs(fitYPoly - origYPoly)
671 # Check that input detector model matches fit
672 np.testing.assert_array_less(absDiffX, 1e-7)
673 np.testing.assert_array_less(absDiffY, 1e-7)
675 def test_missingWcs(self):
676 """Test that task does not fail when the input WCS is None for one
677 extension and that the fit WCS for that extension returns a finite
678 result.
679 """
680 # Need to copy the catalogs since they are modified later.
681 inputVisitSummary = [cat.copy(deep=True) for cat in self.inputVisitSummary]
682 # Set one WCS to be None
683 testVisit = 0
684 testDetector = 20
685 inputVisitSummary[testVisit][testDetector].setWcs(None)
687 outputs = self.task.run(
688 self.inputCatalogRefs,
689 inputVisitSummary,
690 instrumentName=self.instrumentName,
691 refEpoch=self.refEpoch,
692 refObjectLoader=self.refObjectLoader,
693 )
695 # Check that the fit WCS for the extension with input WCS=None returns
696 # finite sky values.
697 testWcs = outputs.outputWcss[self.testVisits[testVisit]][testDetector].getWcs()
698 testSky = testWcs.pixelToSky(0, 0)
699 self.assertTrue(testSky.isFinite())
701 def test_inputCameraModel(self):
702 """Test running task with an input camera model, and check that true
703 object coordinates are recovered.
704 """
705 config = copy(self.config)
706 config.saveCameraModel = False
707 config.useInputCameraModel = True
708 task = GbdesAstrometricFitTask(config=config)
709 with open(os.path.join(self.datadir, "sample_camera_model.yaml"), "r") as f:
710 cameraModel = yaml.load(f, Loader=yaml.Loader)
712 outputs = task.run(
713 self.inputCatalogRefs,
714 self.inputVisitSummary,
715 instrumentName=self.instrumentName,
716 refObjectLoader=self.refObjectLoader,
717 inputCameraModel=cameraModel,
718 )
720 # Check that the output WCS is close (not necessarily exactly equal) to
721 # the result when fitting the full model.
722 for v, visit in enumerate(self.testVisits):
723 visitSummary = self.inputVisitSummary[v]
724 outputWcsCatalog = outputs.outputWcss[visit]
725 visitSources = self.inputCatalogRefs[v].get()
726 for d, detectorRow in enumerate(visitSummary):
727 detectorId = detectorRow["id"]
728 fitwcs = outputWcsCatalog[d].getWcs()
729 detSources = visitSources[visitSources["detector"] == detectorId]
730 fitRA, fitDec = fitwcs.pixelToSkyArray(detSources["x"], detSources["y"], degrees=True)
731 dRA = fitRA - detSources["trueRA"]
732 dDec = fitDec - detSources["trueDec"]
733 # Check that input coordinates match the output coordinates
734 self.assertAlmostEqual(np.mean(dRA), 0)
735 self.assertAlmostEqual(np.std(dRA), 0)
736 self.assertAlmostEqual(np.mean(dDec), 0)
737 self.assertAlmostEqual(np.std(dDec), 0)
739 def test_useColor(self):
740 """Test running task with color catalog and DCR fitting."""
741 config = copy(self.config)
742 config.useColor = True
744 task = GbdesAstrometricFitTask(config=config)
745 outputs = task.run(
746 self.inputCatalogRefs,
747 self.inputVisitSummary,
748 instrumentName=self.instrumentName,
749 refObjectLoader=self.refObjectLoader,
750 colorCatalog=self.colorCatalog,
751 )
752 self.assertEqual(len(outputs.colorParams["visit"]), len(self.testVisits))
755class TestGbdesGlobalAstrometricFit(TestGbdesAstrometricFit):
756 @classmethod
757 def setUpClass(cls):
758 # Set random seed
759 np.random.seed(1234)
761 # Fraction of simulated stars in the reference catalog and science
762 # exposures
763 inReferenceFraction = 1
764 inScienceFraction = 1
766 # Make fake data
767 cls.datadir = os.path.join(TESTDIR, "data")
769 cls.nFields = 2
770 cls.instrumentName = "HSC"
771 cls.refEpoch = 57205.5
773 # Make test inputVisitSummary. VisitSummaryTables are taken from
774 # collection HSC/runs/RC2/w_2022_20/DM-34794
775 cls.testVisits = [1176, 17900, 17930, 17934, 36434, 36460, 36494, 36446]
776 cls.inputVisitSummary = []
777 for testVisit in cls.testVisits:
778 visSum = afwTable.ExposureCatalog.readFits(
779 os.path.join(cls.datadir, f"visitSummary_{testVisit}.fits")
780 )
781 cls.inputVisitSummary.append(visSum)
783 cls.config = GbdesGlobalAstrometricFitConfig()
784 cls.config.systematicError = 0
785 cls.config.devicePolyOrder = 4
786 cls.config.exposurePolyOrder = 6
787 cls.config.fitReserveFraction = 0
788 cls.config.fitReserveRandomSeed = 1234
789 cls.config.saveModelParams = True
790 cls.config.saveCameraModel = True
791 cls.task = GbdesGlobalAstrometricFitTask(config=cls.config)
793 cls.fields, cls.fieldRegions = cls.task._prep_sky(cls.inputVisitSummary)
795 (
796 cls.exposureInfo,
797 cls.exposuresHelper,
798 cls.extensionInfo,
799 cls.instruments,
800 ) = cls.task._get_exposure_info(cls.inputVisitSummary, fieldRegions=cls.fieldRegions)
802 refDataIds, deferredRefCats = [], []
803 allStarIds = []
804 allStarRAs = []
805 allStarDecs = []
806 for region in cls.fieldRegions.values():
807 # Bounding box of observations:
808 bbox = region.getBoundingBox()
809 raMin = bbox.getLon().getA().asDegrees()
810 raMax = bbox.getLon().getB().asDegrees()
811 decMin = bbox.getLat().getA().asDegrees()
812 decMax = bbox.getLat().getB().asDegrees()
814 # Make random set of data in a bounding box determined by input
815 # visits
816 cls.nStars, starIds, starRAs, starDecs = cls._make_simulated_stars(
817 raMin, raMax, decMin, decMax, cls.config.matchRadius
818 )
820 allStarIds.append(starIds)
821 allStarRAs.append(starRAs)
822 allStarDecs.append(starDecs)
824 corners = [
825 lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector()
826 for (ra, dec) in zip([raMin, raMax, raMax, raMin], [decMin, decMin, decMax, decMax])
827 ]
828 conv_box = lsst.sphgeom.ConvexPolygon.convexHull(corners)
829 # Make a reference catalog that will be loaded into
830 # ReferenceObjectLoader
831 refDataId, deferredRefCat = cls._make_refCat(
832 starIds, starRAs, starDecs, inReferenceFraction, conv_box
833 )
834 refDataIds.append(refDataId)
835 deferredRefCats.append(deferredRefCat)
837 cls.refObjectLoader = ReferenceObjectLoader(refDataIds, deferredRefCats)
838 cls.refObjectLoader.config.requireProperMotion = False
839 cls.refObjectLoader.config.anyFilterMapsToThis = "test_filter"
841 cls.task.refObjectLoader = cls.refObjectLoader
843 allRefObjects, allRefCovariances = {}, {}
844 for f, fieldRegion in cls.fieldRegions.items():
845 refObjects, refCovariance = cls.task._load_refcat(
846 cls.refObjectLoader, cls.extensionInfo, epoch=cls.exposureInfo.medianEpoch, region=fieldRegion
847 )
848 allRefObjects[f] = refObjects
849 allRefCovariances[f] = refCovariance
850 cls.refObjects = allRefObjects
852 # Get True WCS for stars:
853 with open(os.path.join(cls.datadir, "sample_global_wcs.yaml"), "r") as f:
854 cls.trueModel = yaml.load(f, Loader=yaml.Loader)
856 cls.trueWCSs = cls._make_wcs(cls.trueModel, cls.inputVisitSummary)
858 cls.isolatedStarCatalogs, cls.isolatedStarSources = cls._make_isolatedStars(
859 allStarIds, allStarRAs, allStarDecs, cls.trueWCSs, inScienceFraction
860 )
862 cls.colorCatalog = cls._make_colors(np.concatenate(allStarRAs), np.concatenate(allStarDecs))
864 cls.outputs = cls.task.run(
865 cls.inputVisitSummary,
866 cls.isolatedStarSources,
867 cls.isolatedStarCatalogs,
868 instrumentName=cls.instrumentName,
869 refEpoch=cls.refEpoch,
870 refObjectLoader=cls.refObjectLoader,
871 )
873 @classmethod
874 def _make_isolatedStars(cls, allStarIds, allStarRAs, allStarDecs, trueWCSs, inScienceFraction):
875 """Given a subset of the simulated data, make source catalogs and star
876 catalogs.
878 This takes the true WCSs to go from the RA and Decs of the simulated
879 stars to pixel coordinates for a given visit and detector. If those
880 pixel coordinates are within the bounding box of the detector, the
881 source and visit information is put in the corresponding catalog of
882 isolated star sources.
884 Parameters
885 ----------
886 allStarIds : `np.ndarray` [`int`]
887 Source ids for the simulated stars
888 allStarRas : `np.ndarray` [`float`]
889 RAs of the simulated stars
890 allStarDecs : `np.ndarray` [`float`]
891 Decs of the simulated stars
892 trueWCSs : `list` [`lsst.afw.geom.SkyWcs`]
893 WCS with which to simulate the source pixel coordinates
894 inReferenceFraction : float
895 Percentage of simulated stars to include in reference catalog
897 Returns
898 -------
899 isolatedStarCatalogRefs : `list`
900 [`lsst.pipe.base.InMemoryDatasetHandle`]
901 List of references to isolated star catalogs.
902 isolatedStarSourceRefs : `list`
903 [`lsst.pipe.base.InMemoryDatasetHandle`]
904 List of references to isolated star sources.
905 """
906 bbox = lsst.geom.BoxD(
907 lsst.geom.Point2D(
908 cls.inputVisitSummary[0][0]["bbox_min_x"], cls.inputVisitSummary[0][0]["bbox_min_y"]
909 ),
910 lsst.geom.Point2D(
911 cls.inputVisitSummary[0][0]["bbox_max_x"], cls.inputVisitSummary[0][0]["bbox_max_y"]
912 ),
913 )
914 bboxCorners = bbox.getCorners()
916 isolatedStarCatalogRefs = []
917 isolatedStarSourceRefs = []
918 for i in range(len(allStarIds)):
919 starIds = allStarIds[i]
920 starRAs = allStarRAs[i]
921 starDecs = allStarDecs[i]
922 isolatedStarCatalog = pd.DataFrame({"ra": starRAs, "dec": starDecs}, index=starIds)
923 isolatedStarCatalogRefs.append(
924 InMemoryDatasetHandle(isolatedStarCatalog, storageClass="DataFrame", dataId={"tract": 0})
925 )
926 sourceCats = []
927 for v, visit in enumerate(cls.testVisits):
928 nVisStars = int(cls.nStars * inScienceFraction)
929 visitStarIndices = np.random.choice(cls.nStars, nVisStars, replace=False)
930 visitStarIds = starIds[visitStarIndices]
931 visitStarRas = starRAs[visitStarIndices]
932 visitStarDecs = starDecs[visitStarIndices]
933 for detector in trueWCSs[visit]:
934 detWcs = detector.getWcs()
935 detectorId = detector["id"]
936 radecCorners = detWcs.pixelToSky(bboxCorners)
937 detectorFootprint = sphgeom.ConvexPolygon([rd.getVector() for rd in radecCorners])
938 detectorIndices = detectorFootprint.contains(
939 (visitStarRas * u.degree).to(u.radian), (visitStarDecs * u.degree).to(u.radian)
940 )
941 nDetectorStars = detectorIndices.sum()
942 detectorArray = np.ones(nDetectorStars, dtype=int) * detector["id"]
943 visitArray = np.ones(nDetectorStars, dtype=int) * visit
945 ones_like = np.ones(nDetectorStars)
946 zeros_like = np.zeros(nDetectorStars, dtype=bool)
948 x, y = detWcs.skyToPixelArray(
949 visitStarRas[detectorIndices], visitStarDecs[detectorIndices], degrees=True
950 )
952 origWcs = (cls.inputVisitSummary[v][cls.inputVisitSummary[v]["id"] == detectorId])[
953 0
954 ].getWcs()
955 inputRa, inputDec = origWcs.pixelToSkyArray(x, y, degrees=True)
957 sourceDict = {}
958 sourceDict["detector"] = detectorArray
959 sourceDict["visit"] = visitArray
960 sourceDict["obj_index"] = visitStarIds[detectorIndices]
961 sourceDict["x"] = x
962 sourceDict["y"] = y
963 sourceDict["xErr"] = 1e-3 * ones_like
964 sourceDict["yErr"] = 1e-3 * ones_like
965 sourceDict["inputRA"] = inputRa
966 sourceDict["inputDec"] = inputDec
967 sourceDict["trueRA"] = visitStarRas[detectorIndices]
968 sourceDict["trueDec"] = visitStarDecs[detectorIndices]
969 for key in ["apFlux_12_0_flux", "apFlux_12_0_instFlux", "ixx", "iyy"]:
970 sourceDict[key] = ones_like
971 for key in [
972 "pixelFlags_edge",
973 "pixelFlags_saturated",
974 "pixelFlags_interpolatedCenter",
975 "pixelFlags_interpolated",
976 "pixelFlags_crCenter",
977 "pixelFlags_bad",
978 "hsmPsfMoments_flag",
979 "apFlux_12_0_flag",
980 "extendedness",
981 "parentSourceId",
982 "deblend_nChild",
983 "ixy",
984 ]:
985 sourceDict[key] = zeros_like
986 sourceDict["apFlux_12_0_instFluxErr"] = 1e-2 * ones_like
987 sourceDict["detect_isPrimary"] = ones_like.astype(bool)
989 sourceCat = pd.DataFrame(sourceDict)
990 sourceCats.append(sourceCat)
992 isolatedStarSourceTable = pd.concat(sourceCats, ignore_index=True)
993 isolatedStarSourceTable = isolatedStarSourceTable.sort_values(by=["obj_index"])
994 isolatedStarSourceRefs.append(
995 InMemoryDatasetHandle(isolatedStarSourceTable, storageClass="DataFrame", dataId={"tract": 0})
996 )
998 return isolatedStarCatalogRefs, isolatedStarSourceRefs
1000 def test_loading_and_association(self):
1001 """Test that associated objects actually correspond to the same
1002 simulated object."""
1003 associations, sourceDict = self.task._associate_from_isolated_sources(
1004 self.isolatedStarSources, self.isolatedStarCatalogs, self.extensionInfo, self.refObjects
1005 )
1007 object_groups = np.flatnonzero(np.array(associations.sequence) == 0)
1008 for i in range(len(object_groups) - 1)[:10]:
1009 ras, decs = [], []
1010 for ind in np.arange(object_groups[i], object_groups[i + 1]):
1011 visit = self.extensionInfo.visit[associations.extn[ind]]
1012 detectorInd = self.extensionInfo.detectorIndex[associations.extn[ind]]
1013 detector = self.extensionInfo.detector[associations.extn[ind]]
1014 if detectorInd == -1:
1015 ra = self.refObjects[visit * -1]["ra"][associations.obj[ind]]
1016 dec = self.refObjects[visit * -1]["dec"][associations.obj[ind]]
1017 ras.append(ra)
1018 decs.append(dec)
1019 else:
1020 x = sourceDict[visit][detector]["x"][associations.obj[ind]]
1021 y = sourceDict[visit][detector]["y"][associations.obj[ind]]
1022 ra, dec = self.trueWCSs[visit][detectorInd].getWcs().pixelToSky(x, y)
1023 ras.append(ra.asDegrees())
1024 decs.append(dec.asDegrees())
1025 np.testing.assert_allclose(ras, ras[0])
1026 np.testing.assert_allclose(decs, decs[0])
1028 def test_refCatLoader(self):
1029 """Test loading objects from the refCat in each of the fields."""
1031 for region in self.fieldRegions.values():
1032 refCat, refCov = self.task._load_refcat(
1033 self.refObjectLoader, self.extensionInfo, region=region, epoch=self.exposureInfo.medianEpoch
1034 )
1035 assert len(refCat) > 0
1037 def test_make_outputs(self):
1038 """Test that the run method recovers the input model parameters."""
1039 for isolatedStarSourceRef in self.isolatedStarSources:
1040 iss = isolatedStarSourceRef.get()
1041 visits = np.unique(iss["visit"])
1042 for v, visit in enumerate(visits):
1043 outputWcsCatalog = self.outputs.outputWcss[visit]
1044 visitSources = iss[iss["visit"] == visit]
1045 detectors = outputWcsCatalog["id"]
1046 for d, detectorId in enumerate(detectors):
1047 fitwcs = outputWcsCatalog[d].getWcs()
1048 detSources = visitSources[visitSources["detector"] == detectorId]
1049 fitRA, fitDec = fitwcs.pixelToSkyArray(detSources["x"], detSources["y"], degrees=True)
1050 dRA = fitRA - detSources["trueRA"]
1051 dDec = fitDec - detSources["trueDec"]
1052 # Check that input coordinates match the output coordinates
1053 self.assertAlmostEqual(np.mean(dRA), 0)
1054 self.assertAlmostEqual(np.std(dRA), 0)
1055 self.assertAlmostEqual(np.mean(dDec), 0)
1056 self.assertAlmostEqual(np.std(dDec), 0)
1058 def test_missingWcs(self):
1059 """Test that task does not fail when the input WCS is None for one
1060 extension and that the fit WCS for that extension returns a finite
1061 result.
1062 """
1063 # Need to copy the catalogs since they are modified later.
1064 inputVisitSummary = [cat.copy(deep=True) for cat in self.inputVisitSummary]
1065 # Set one WCS to be None
1066 testVisit = 0
1067 testDetector = 20
1068 inputVisitSummary[testVisit][testDetector].setWcs(None)
1070 outputs = self.task.run(
1071 inputVisitSummary,
1072 self.isolatedStarSources,
1073 self.isolatedStarCatalogs,
1074 instrumentName=self.instrumentName,
1075 refEpoch=self.refEpoch,
1076 refObjectLoader=self.refObjectLoader,
1077 )
1079 # Check that the fit WCS for the extension with input WCS=None returns
1080 # finite sky values.
1081 testWcs = outputs.outputWcss[self.testVisits[testVisit]][testDetector].getWcs()
1082 testSky = testWcs.pixelToSky(0, 0)
1083 self.assertTrue(testSky.isFinite())
1085 def test_inputCameraModel(self):
1086 """Test running task with an input camera model, and check that true
1087 object coordinates are recovered.
1088 """
1089 config = copy(self.config)
1090 config.saveCameraModel = False
1091 config.useInputCameraModel = True
1092 task = GbdesGlobalAstrometricFitTask(config=config)
1093 with open(os.path.join(self.datadir, "sample_global_camera_model.yaml"), "r") as f:
1094 cameraModel = yaml.load(f, Loader=yaml.Loader)
1096 outputs = task.run(
1097 self.inputVisitSummary,
1098 self.isolatedStarSources,
1099 self.isolatedStarCatalogs,
1100 instrumentName=self.instrumentName,
1101 refObjectLoader=self.refObjectLoader,
1102 inputCameraModel=cameraModel,
1103 )
1105 # Check that the output WCS is close (not necessarily exactly equal) to
1106 # the result when fitting the full model.
1107 for isolatedStarSourceRef in self.isolatedStarSources:
1108 iss = isolatedStarSourceRef.get()
1109 visits = np.unique(iss["visit"])
1110 for v, visit in enumerate(visits):
1111 outputWcsCatalog = outputs.outputWcss[visit]
1112 visitSources = iss[iss["visit"] == visit]
1113 detectors = outputWcsCatalog["id"]
1114 for d, detectorId in enumerate(detectors):
1115 fitwcs = outputWcsCatalog[d].getWcs()
1116 detSources = visitSources[visitSources["detector"] == detectorId]
1117 fitRA, fitDec = fitwcs.pixelToSkyArray(detSources["x"], detSources["y"], degrees=True)
1118 dRA = fitRA - detSources["trueRA"]
1119 dDec = fitDec - detSources["trueDec"]
1120 # Check that input coordinates match the output coordinates
1121 self.assertAlmostEqual(np.mean(dRA), 0, places=6)
1122 self.assertAlmostEqual(np.std(dRA), 0)
1123 self.assertAlmostEqual(np.mean(dDec), 0, places=6)
1124 self.assertAlmostEqual(np.std(dDec), 0)
1126 def test_useColor(self):
1127 """Test running task with color catalog and DCR fitting."""
1128 config = copy(self.config)
1129 config.useColor = True
1131 task = GbdesGlobalAstrometricFitTask(config=config)
1133 outputs = task.run(
1134 self.inputVisitSummary,
1135 self.isolatedStarSources,
1136 self.isolatedStarCatalogs,
1137 instrumentName=self.instrumentName,
1138 refObjectLoader=self.refObjectLoader,
1139 colorCatalog=self.colorCatalog,
1140 )
1141 self.assertEqual(len(outputs.colorParams["visit"]), len(self.testVisits))
1144class TestApparentMotion(lsst.utils.tests.TestCase):
1145 """This class simulates an array of objects with random proper motions and
1146 parallaxes and compares the results of
1147 `lsst.drp.tasks.gbdesAstrometricFit.calculate_apparent_motion` against the
1148 results calculated by astropy."""
1150 def setUp(self):
1151 np.random.seed(12345)
1153 self.ras = (np.random.rand(10) + 150.0) * u.degree
1154 self.decs = (np.random.rand(10) + 2.5) * u.degree
1156 self.pmRAs = (np.random.random(10) * 10) * u.mas / u.yr
1157 self.pmDecs = (np.random.random(10) * 10) * u.mas / u.yr
1159 self.parallaxes = (abs(np.random.random(10) * 5)) * u.mas
1161 self.mjds = astropy.time.Time(np.random.rand(10) * 20 + 57000, format="mjd", scale="tai")
1163 self.refEpoch = astropy.time.Time(57100, format="mjd", scale="tai")
1165 def test_proper_motion(self):
1166 """Calculate the change in position due to proper motion only for a
1167 given time change, and compare the results against astropy.
1168 """
1169 # Set parallaxes to zero to test proper motion part by itself.
1170 parallaxes = np.zeros(len(self.ras)) * u.mas
1172 data = astropy.table.Table(
1173 {
1174 "ra": self.ras,
1175 "dec": self.decs,
1176 "pmRA": self.pmRAs,
1177 "pmDec": self.pmDecs,
1178 "parallax": parallaxes,
1179 "MJD": self.mjds,
1180 }
1181 )
1182 raCorrection, decCorrection = calculate_apparent_motion(data, self.refEpoch)
1184 pmRACosDec = self.pmRAs * np.cos(self.decs)
1185 coords = SkyCoord(
1186 ra=self.ras,
1187 dec=self.decs,
1188 pm_ra_cosdec=pmRACosDec,
1189 pm_dec=self.pmDecs,
1190 obstime=self.refEpoch,
1191 )
1193 newCoords = coords.apply_space_motion(new_obstime=self.mjds)
1194 astropyDRA = newCoords.ra.degree - coords.ra.degree
1195 astropyDDec = newCoords.dec.degree - coords.dec.degree
1197 # The proper motion values are on the order of 0-3 mas/yr and the
1198 # differences are on the scale of 1e-8 mas.
1199 self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=1e-6)
1200 self.assertFloatsAlmostEqual(astropyDDec, decCorrection.value, rtol=1e-6)
1202 def test_parallax(self):
1203 """Calculate the change in position due to parallax for a given time
1204 change and compare the results against astropy. Astropy will not give
1205 the parallax part separate from other sources of space motion, such as
1206 annual aberration, but it can be obtained by comparing the calculated
1207 positions of an object with a given parallax and one with a much
1208 further distance. The result is still only approximately the same, but
1209 is close enough for accuracy needed here.
1210 """
1211 # Set proper motions to zero to test parallax correction by itself.
1212 pms = np.zeros(len(self.ras)) * u.mas / u.yr
1214 data = astropy.table.Table(
1215 {
1216 "ra": self.ras,
1217 "dec": self.decs,
1218 "pmRA": pms,
1219 "pmDec": pms,
1220 "parallax": self.parallaxes,
1221 "MJD": self.mjds,
1222 }
1223 )
1224 raCorrection, decCorrection = calculate_apparent_motion(data, self.refEpoch)
1226 coords = SkyCoord(
1227 ra=self.ras,
1228 dec=self.decs,
1229 pm_ra_cosdec=pms,
1230 pm_dec=pms,
1231 distance=Distance(parallax=self.parallaxes),
1232 obstime=self.refEpoch,
1233 )
1234 # Make another set of objects at the same coordinates but a much
1235 # greater distant in order to isolate parallax when applying space
1236 # motion.
1237 refCoords = SkyCoord(
1238 ra=self.ras,
1239 dec=self.decs,
1240 pm_ra_cosdec=pms,
1241 pm_dec=pms,
1242 distance=1e10 * Distance(parallax=self.parallaxes),
1243 obstime=self.refEpoch,
1244 )
1246 newCoords = coords.apply_space_motion(new_obstime=self.mjds).transform_to("gcrs")
1247 refNewCoords = refCoords.apply_space_motion(new_obstime=self.mjds).transform_to("gcrs")
1248 astropyDRA = newCoords.ra.degree - refNewCoords.ra.degree
1249 astropyDDec = newCoords.dec.degree - refNewCoords.dec.degree
1251 # The parallax values are on the scale of 1-3 mas, and the difference
1252 # between the two calculation methods is ~0.05 mas.
1253 self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=2e-2)
1254 self.assertFloatsAlmostEqual(astropyDDec, decCorrection.value, rtol=2e-2)
1257def setup_module(module):
1258 lsst.utils.tests.init()
1261if __name__ == "__main__": 1261 ↛ 1262line 1261 didn't jump to line 1262 because the condition on line 1261 was never true
1262 lsst.utils.tests.init()
1263 unittest.main()