882def applyProperMotionsImpl(log, catalog, epoch):
883 """Apply proper motion correction to a reference catalog.
884
885 Adjust position and position error in the ``catalog``
886 for proper motion to the specified ``epoch``,
887 modifying the catalog in place.
888
889 Parameters
890 ----------
891 log : `lsst.log.Log` or `logging.Logger`
892 Log object to write to.
893 catalog : `lsst.afw.table.SimpleCatalog`
894 Catalog of positions, containing:
895
896 - Coordinates, retrieved by the table's coordinate key.
897 - ``coord_raErr`` : Error in Right Ascension (rad).
898 - ``coord_decErr`` : Error in Declination (rad).
899 - ``pm_ra`` : Proper motion in Right Ascension (rad/yr,
900 East positive)
901 - ``pm_raErr`` : Error in ``pm_ra`` (rad/yr), optional.
902 - ``pm_dec`` : Proper motion in Declination (rad/yr,
903 North positive)
904 - ``pm_decErr`` : Error in ``pm_dec`` (rad/yr), optional.
905 - ``epoch`` : Mean epoch of object (an astropy.time.Time)
906 epoch : `astropy.time.Time`
907 Epoch to which to correct proper motion.
908 """
909 if "epoch" not in catalog.schema or "pm_ra" not in catalog.schema or "pm_dec" not in catalog.schema:
910 log.warning("Proper motion correction not available from catalog")
911 return
912 if not catalog.isContiguous():
913 raise RuntimeError("Catalog must be contiguous")
914 catEpoch = astropy.time.Time(catalog["epoch"], scale="tai", format="mjd")
915 log.info("Correcting reference catalog for proper motion to %r", epoch)
916
917 timeDiffsYears = (epoch.tai - catEpoch).to(astropy.units.yr).value
918 coordKey = catalog.table.getCoordKey()
919
920
921 pmRaRad = catalog["pm_ra"]
922 pmDecRad = catalog["pm_dec"]
923 offsetsRaRad = pmRaRad*timeDiffsYears
924 offsetsDecRad = pmDecRad*timeDiffsYears
925
926
927
928
929
930
931
932 offsetBearingsRad = numpy.arctan2(offsetsDecRad*1e6, offsetsRaRad*1e6)
933 offsetAmountsRad = numpy.hypot(offsetsRaRad, offsetsDecRad)
934 for record, bearingRad, amountRad in zip(catalog, offsetBearingsRad, offsetAmountsRad):
935 record.set(coordKey,
936 record.get(coordKey).offset(bearing=bearingRad*geom.radians,
937 amount=amountRad*geom.radians))
938
939
940 if "coord_raErr" in catalog.schema:
941 catalog["coord_raErr"] = numpy.hypot(catalog["coord_raErr"],
942 catalog["pm_raErr"]*timeDiffsYears)
943 if "coord_decErr" in catalog.schema:
944 catalog["coord_decErr"] = numpy.hypot(catalog["coord_decErr"],
945 catalog["pm_decErr"]*timeDiffsYears)