833def applyProperMotionsImpl(log, catalog, epoch):
834 """Apply proper motion correction to a reference catalog.
835
836 Adjust position and position error in the ``catalog``
837 for proper motion to the specified ``epoch``,
838 modifying the catalog in place.
839
840 Parameters
841 ----------
842 log : `lsst.log.Log` or `logging.getLogger`
843 Log object to write to.
844 catalog : `lsst.afw.table.SimpleCatalog`
845 Catalog of positions, containing:
846
847 - Coordinates, retrieved by the table's coordinate key.
848 - ``coord_raErr`` : Error in Right Ascension (rad).
849 - ``coord_decErr`` : Error in Declination (rad).
850 - ``pm_ra`` : Proper motion in Right Ascension (rad/yr,
851 East positive)
852 - ``pm_raErr`` : Error in ``pm_ra`` (rad/yr), optional.
853 - ``pm_dec`` : Proper motion in Declination (rad/yr,
854 North positive)
855 - ``pm_decErr`` : Error in ``pm_dec`` (rad/yr), optional.
856 - ``epoch`` : Mean epoch of object (an astropy.time.Time)
857 epoch : `astropy.time.Time`
858 Epoch to which to correct proper motion.
859 """
860 if "epoch" not in catalog.schema or "pm_ra" not in catalog.schema or "pm_dec" not in catalog.schema:
861 log.warning("Proper motion correction not available from catalog")
862 return
863 if not catalog.isContiguous():
864 raise RuntimeError("Catalog must be contiguous")
865 catEpoch = astropy.time.Time(catalog["epoch"], scale="tai", format="mjd")
866 log.info("Correcting reference catalog for proper motion to %r", epoch)
867
868 timeDiffsYears = (epoch.tai - catEpoch).to(astropy.units.yr).value
869 coordKey = catalog.table.getCoordKey()
870
871
872 pmRaRad = catalog["pm_ra"]
873 pmDecRad = catalog["pm_dec"]
874 offsetsRaRad = pmRaRad*timeDiffsYears
875 offsetsDecRad = pmDecRad*timeDiffsYears
876
877
878
879
880
881
882
883 offsetBearingsRad = numpy.arctan2(offsetsDecRad*1e6, offsetsRaRad*1e6)
884 offsetAmountsRad = numpy.hypot(offsetsRaRad, offsetsDecRad)
885 for record, bearingRad, amountRad in zip(catalog, offsetBearingsRad, offsetAmountsRad):
886 record.set(coordKey,
887 record.get(coordKey).offset(bearing=bearingRad*geom.radians,
888 amount=amountRad*geom.radians))
889
890
891 if "coord_raErr" in catalog.schema:
892 catalog["coord_raErr"] = numpy.hypot(catalog["coord_raErr"],
893 catalog["pm_raErr"]*timeDiffsYears)
894 if "coord_decErr" in catalog.schema:
895 catalog["coord_decErr"] = numpy.hypot(catalog["coord_decErr"],
896 catalog["pm_decErr"]*timeDiffsYears)