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