Source code for lsst.sims.maf.stackers.coordStackers

import numpy as np
import ephem
from lsst.sims.utils import _galacticFromEquatorial, calcLmstLast

from .baseStacker import BaseStacker
from .ditherStackers import wrapRA

__all__ = ['mjd2djd', 'raDec2AltAz', 'GalacticStacker', 'EclipticStacker']


[docs]def mjd2djd(mjd): """Convert MJD to the Dublin Julian date used by ephem. Parameters ---------- mjd : float or numpy.ndarray The modified julian date. Returns ------- float or numpy.ndarray The dublin julian date. """ doff = ephem.Date(0)-ephem.Date('1858/11/17') djd = mjd-doff return djd
[docs]def raDec2AltAz(ra, dec, lat, lon, mjd, altonly=False): """Convert RA/Dec (and telescope site lat/lon) to alt/az. This uses simple equations and ignores aberation, precession, nutation, etc. Parameters ---------- ra : array_like RA, in radians. dec : array_like Dec, in radians. Must be same length as `ra`. lat : float Latitude of the observatory in radians. lon : float Longitude of the observatory in radians. mjd : float Modified Julian Date. altonly : bool, opt Calculate altitude only. Returns ------- alt : numpy.array Altitude, same length as `ra` and `dec`. Radians. az : numpy.array Azimuth, same length as `ra` and `dec`. Radians. """ lmst, last = calcLmstLast(mjd, lon) lmst = lmst / 12. * np.pi # convert to rad ha = lmst - ra sindec = np.sin(dec) sinlat = np.sin(lat) coslat = np.cos(lat) sinalt = sindec * sinlat + np.cos(dec) * coslat * np.cos(ha) # make sure sinalt is in the expected range. sinalt = np.where(sinalt < -1, -1, sinalt) sinalt = np.where(sinalt > 1, 1, sinalt) alt = np.arcsin(sinalt) if altonly: az = None else: cosaz = (sindec-np.sin(alt)*sinlat)/(np.cos(alt)*coslat) cosaz = np.where(cosaz < -1, -1, cosaz) cosaz = np.where(cosaz > 1, 1, cosaz) az = np.arccos(cosaz) signflip = np.where(np.sin(ha) > 0) az[signflip] = 2.*np.pi-az[signflip] return alt, az
[docs]class GalacticStacker(BaseStacker): """Add the galactic coordinates of each RA/Dec pointing: gall, galb Parameters ---------- raCol : str, opt Name of the RA column. Default fieldRA. decCol : str, opt Name of the Dec column. Default fieldDec. """ colsAdded = ['gall', 'galb'] def __init__(self, raCol='fieldRA', decCol='fieldDec', degrees=True): self.colsReq = [raCol, decCol] self.raCol = raCol self.decCol = decCol self.degrees = degrees if self.degrees: self.units = ['degrees', 'degrees'] else: self.units = ['radians', 'radians'] def _run(self, simData, cols_present=False): # raCol and DecCol in radians, gall/b in radians. if cols_present: # Column already present in data; assume it is correct and does not need recalculating. return simData if self.degrees: simData['gall'], simData['galb'] = _galacticFromEquatorial(np.radians(simData[self.raCol]), np.radians(simData[self.decCol])) else: simData['gall'], simData['galb'] = _galacticFromEquatorial(simData[self.raCol], simData[self.decCol]) return simData
[docs]class EclipticStacker(BaseStacker): """Add the ecliptic coordinates of each RA/Dec pointing: eclipLat, eclipLon Optionally subtract off the sun's ecliptic longitude and wrap. Parameters ---------- mjdCol : str, opt Name of the MJD column. Default expMJD. raCol : str, opt Name of the RA column. Default fieldRA. decCol : str, opt Name of the Dec column. Default fieldDec. subtractSunLon : bool, opt Flag to subtract the sun's ecliptic longitude. Default False. """ colsAdded = ['eclipLat', 'eclipLon'] def __init__(self, mjdCol='observationStartMJD', raCol='fieldRA', decCol='fieldDec', degrees=True, subtractSunLon=False): self.colsReq = [mjdCol, raCol, decCol] self.subtractSunLon = subtractSunLon self.degrees = degrees if self.degrees: self.units = ['degrees', 'degrees'] else: self.units = ['radians', 'radians'] self.mjdCol = mjdCol self.raCol = raCol self.decCol = decCol def _run(self, simData, cols_present=False): if cols_present: # Column already present in data; assume it is correct and does not need recalculating. return simData for i in np.arange(simData.size): if self.degrees: coord = ephem.Equatorial(np.radians(simData[self.raCol][i]), np.radians(simData[self.decCol][i]), epoch=2000) else: coord = ephem.Equatorial(simData[self.raCol][i], simData[self.decCol][i], epoch=2000) ecl = ephem.Ecliptic(coord) simData['eclipLat'][i] = ecl.lat if self.subtractSunLon: djd = mjd2djd(simData[self.mjdCol][i]) sun = ephem.Sun(djd) sunEcl = ephem.Ecliptic(sun) lon = wrapRA(ecl.lon - sunEcl.lon) simData['eclipLon'][i] = lon else: simData['eclipLon'][i] = ecl.lon if self.degrees: simData['eclipLon'] = np.degrees(simData['eclipLon']) simData['eclipLat'] = np.degrees(simData['eclipLat']) return simData