# Looking for a routine to correct for Local Standard of Rest (LSR)

I've been scanning the Milky Way with a Small Single Dish Radio Telescope, obtaining Spectral data with a Span of ~ 2 MHz, centered at the 21 cm line. With this information I can derive a relative velocity (due to Doppler shifted emission) between the observer (telescope, at lat ~ 33° South) and gas clouds in the line of sight. To study the position of this clouds relative to the Galactic Center, I know that the way to correct for this is subtracting our motion relative to this 'fixed' point: correcting by Earth's rotation, Earth motion on the Solar System plane and the System's motion itself (everything projected onto the line of sight), so my data is comparable to any other observation.

I've been reading on related questions such as

I also looked on "Tools of Radio Astronomy", T.L Wilson, K. Rohlfs, S. Hüttemeister, 5ed., pg. 189. that briefly describes this procedure.  So I wonder how this is really done in practice and if there's any, e.g., Python Package that does it (maybe in Astropy?).

EDIT: Thanks to @ELNJ 's response and help I managed to get agreeable results. This is data from a Standard Region, particularly S9.

def Freq_to_v(freqval, freq0, mjd_time = 59244.8, obs_ra = 268.1, obs_dec = -34.43):

#This is the observation for S9, 30-01-2021, 16:34 hrs CL.
#S9 Std. Reg. Obs. time = 59244.8

t = Time(mjd_time, format='mjd', scale = 'utc')

#Obs. Location and height in GEODETIC coordinates (and not geocentric a X Y Z quantities tuple!!).

loc = EarthLocation(lon= 289.47, lat =-33.27, height = 1450)

sc = SkyCoord(1*u.deg, 2*u.deg)

#####################################################

#Regular Doppler Shift Calculation. Here are two conventions: optical and radio doppler shift.

shift_opt, shift_rad = (freqval - freq0)/freqval, (freq0-freqval)/freq0 #shifting factor, no dimensions.

#################################################### To barycentric:

rv = v_rad - vcorr #corrected to barycentric w/o Special Relativity terms.

#This is the observation for S9 Std, Region., 30-01-2021, 16:34 hrs CL.

################################################### To LSR:

my_observation = ICRS(ra=obs_ra *u.deg, dec= obs_dec *u.deg, \
pm_ra_cosdec=0*u.mas/u.yr, pm_dec=0*u.mas/u.yr, \

return new_rv



You can compare Williams' 1972 result with mine: and check that it's in fact consistent.

There are several corrections here, and as you suspect, all of them can be done in astropy. You'll need to do something like this:

1. Transform your observed velocity to the solar system barycenter (thus correcting for Earth's rotation on its axis, and its orbit around the Sun. There is an example here in the astropy docs that shows how to do that. If you follow that example, you'll first need to define your radial velocity as an astropy "quantity" object with units, e.g. something like this:
import astropy.units as u
# Define RV of 30 km/s
rv = 30 * u.km / u.s


You can ignore the discussion in that example of the relativistic correction, and just add the barycentric correction to your velocity, e.g. corrected_rv = rv + vcorr

1. Then, you can correct for the fact that the Sun isn't moving quite in step with the average of the nearby stars, and convert from the barycentric radial velocity to the local standard of rest, or LSR. Astropy has some examples here but they don't do anything quite like this, so here's some example code.

Set up an observation that has the coordinates you observed, and the corrected radial velocity from step 1. (I'm continuing with 30 km/s here, but you would substitute whatever you got in step 1.)

from astropy.coordinates import ICRS, LSR
my_observation = ICRS(ra=17.75*15*u.deg, dec=-29*u.deg, \
pm_ra_cosdec=0*u.mas/u.yr, pm_dec=0*u.mas/u.yr, \

new_rv = my_observation.transform_to(LSR()).radial_velocity