LAMOST DR5v3

# Clean-up of LAMOST DR5v3 catalog


## LAMOST survey 
LAMOST survey is a China-led stellar spectroscopic survey. The latest catalog, LAMOST DR5 catalog (DR5v3), provides one of the largest spectroscopic catalogs of halo stars ( http://dr5.lamost.org ). 

It contains `5,348,712` stars in total (http://dr5.lamost.org/doc/data-production-description#toc_16). Even if we select a subset of stars with [Fe/H] < -1, we still have `151,954` stars. 

## Radial velocity correction

Due to the systematic error in LAMOST data, we need an empirical correction given by:

- data['rv'] +  4.5  

as the radial velocity. (data['rv'] is the raw data of the radial velocity.)


## Radial velocity random error `rv_err` 

174 stars in DR5 have  `rv_err = -9999`.

Interestingly, the distribution of `rv` of those stars with `rv_err = -9999` does not look peculiar. 

[Strategy] 

- We can discard those stars with `rv_err < 0`. This does not reduce the sample size significantly (174 stars out of 151,954 stars). If our method could find a certian cluster with high confidence, the signal of the cluster should not be dependent on whether we include data with poor `rv` or not. 


![Screen Shot 2020-05-20 at 4.31.02 PM.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610740/844e60d8-1354-ef80-fb39-121f36260ccf.png)


## Duplicated sources 

The stellar parameters (such as [Fe/H] and rv) are derived for each spectrum. There are multiple sources for which the stellar spectra were taken more than once. Such a star appears in the catalog more than once. 

The presence of duplicated sources is problematic when we perform cluster findings because such a star can create a fake signal of a cluster. 

[Strategy] 

- After removing all the entries with `rv_err<0`, we sort the table according to `rv_err`. Then use `astropy` functionality named `unique` to select only the first entry for any given `designation` (target's specfic name; such as `J000001.45+180556.7`) 

- ~~We first crossmatch the LAMOST DR5 against Gaia DR2. If we find multiple entries for a given Gaia DR2 `source_id` (i.e., multiple entries for a unique star), we select an entry with the smallest `rv_err` (assuming that we have discarded stars with `rv_err<0`).~~


---

# The procedure that I followed.

- (1) Download LAMOST DR5v3. 

- (2) Select entries with `feh < -1` and `0 < rv_err < 20`.

- (3) For each `designation` in LAMOST DR5, select an entry with the lowest `rv_err`. 

- (4) Crossmatch against Gaia DR2 by using ADQL script on TOPCAT (cone radius of 1/3600 degree = 1 arcsec). 


```python

#
# ------------
# system package
#
import sys
import os
#
# ------------
# basic package
#
import numpy
import scipy

from astropy import table
from astropy.table import Table
from astropy.table import Column
from astropy.io import fits
from astropy.table import unique
#
# ------------
# astronomical package (PyGaia)
#
from pygaia.astrometry.coordinates import CoordinateTransformation, Transformations
from pygaia.astrometry.vectorastrometry import cartesianToSpherical, astrometryToPhaseSpace, phaseSpaceToAstrometry
from pygaia.errors.astrometric import properMotionErrorSkyAvg, properMotionMinError, properMotionMaxError
from pygaia.photometry.utils import vminiFromSpt, vabsFromSpt
from pygaia.photometry.transformations import gminvFromVmini
# ------------
# astronomical package (Kohei Hattori) -- not cleaned yet.
#
from useful_functions_for_conversion import *
#

#----------

#input_filename = 'dr5v3_stellar.fits'    # This is the data that one needs to use as the first step.
input_filename = 'dr5_stellar_FeHm1.fits' # This is a cheat. I already selected stars with feh<-1.

input_file_full_path = '/Users/kohei/rprocess/NearbyGaiaData/%s' % (input_filename)
hdul = fits.open(input_file_full_path)

#header_list = hdul[0].header
# The header includes some information on who made this file etc. Not important scientifically.

columns = hdul[1].columns
data0   = hdul[1].data

#print(columns.info() )
#print('-'*20)
#print(columns.info('name'))
# >>>
#name:
#    ['obsid', 'designation', 'obsdate', 'lmjd', 'mjd', 'planid', 'spid', 'fiberid', 'ra_obs', 'dec_obs', 'snru', 'snrg', 'snrr', 'snri', 'snrz', 'objtype', 'class', 'subclass', 'z', 'z_err', 'magtype', 'mag1', 'mag2', 'mag3', 'mag4', 'mag5', 'mag6', 'mag7', 'tsource', 'fibertype', 'tfrom', 'tcomment', 'offsets', 'offset_v', 'ra', 'dec', 'teff', 'teff_err', 'logg', 'logg_err', 'feh', 'feh_err', 'rv', 'rv_err']
#None
#print('-'*20)
print(columns.names)
# >>>
#['obsid', 'designation', 'obsdate', 'lmjd', 'mjd', 'planid', 'spid', 'fiberid', 'ra_obs', 'dec_obs', 'snru', 'snrg', 'snrr', 'snri', 'snrz', 'objtype', 'class', 'subclass', 'z', 'z_err', 'magtype', 'mag1', 'mag2', 'mag3', 'mag4', 'mag5', 'mag6', 'mag7', 'tsource', 'fibertype', 'tfrom', 'tcomment', 'offsets', 'offset_v', 'ra', 'dec', 'teff', 'teff_err', 'logg', 'logg_err', 'feh', 'feh_err', 'rv', 'rv_err']

print('-'*20)

#sys.exit(0)

tmp_feh     = data0['feh']
tmp_rv_err  = data0['rv_err']
tmp_ra_deg  = data0['ra'] #['ra_epoch2000'] #['ra']
tmp_dec_deg = data0['dec']#['dec_epoch2000']#['dec']


tmp_mag1 = data0['mag1']
tmp_mag2 = data0['mag2']
tmp_mag3 = data0['mag3']
tmp_mag4 = data0['mag4']
tmp_mag5 = data0['mag5']
tmp_mag6 = data0['mag6']
tmp_mag7 = data0['mag7']


#################################################################
# The photometric information in LAMOST DR5v3 is weird.
# So I just compute the mean magnitude of passbands shown
# (but negleting those magniutdes outside 5-25 mag range).
#
sel_reasonable_m1 = (5<tmp_mag1)*(tmp_mag1<25)
sel_reasonable_m2 = (5<tmp_mag2)*(tmp_mag2<25)
sel_reasonable_m3 = (5<tmp_mag3)*(tmp_mag3<25)
sel_reasonable_m4 = (5<tmp_mag4)*(tmp_mag4<25)
sel_reasonable_m5 = (5<tmp_mag5)*(tmp_mag5<25)
sel_reasonable_m6 = (5<tmp_mag6)*(tmp_mag6<25)
sel_reasonable_m7 = (5<tmp_mag7)*(tmp_mag7<25)

N = len(tmp_mag1)
w1 = numpy.zeros(N) #Initial weight=0
w2 = numpy.zeros(N)
w3 = numpy.zeros(N)
w4 = numpy.zeros(N)
w5 = numpy.zeros(N)
w6 = numpy.zeros(N)
w7 = numpy.zeros(N)

w1[sel_reasonable_m1] = 1 #If reasonable photometric value is given, set weight=1
w2[sel_reasonable_m2] = 1
w3[sel_reasonable_m3] = 1
w4[sel_reasonable_m4] = 1
w5[sel_reasonable_m5] = 1
w6[sel_reasonable_m6] = 1
w7[sel_reasonable_m7] = 1

w_all = (w1+w2+w3+w4+w5+w6+w7) + 1e-6 # number of effective passbands used. (Avoid division by zero.)

weighted_mean = (\
w1*tmp_mag1 + \
w2*tmp_mag2 + \
w3*tmp_mag3 + \
w4*tmp_mag4 + \
w5*tmp_mag5 + \
w6*tmp_mag6 + \
w7*tmp_mag7 ) \
/ w_all

#################################################################
# LAMOST does not provide (l,b). So we create new columns.
# We do not need (pml,pmb), so we use dummy variables named "_".
#
tmp_ra_radian  = deg2radian(tmp_ra_deg)
tmp_dec_radian = deg2radian(tmp_dec_deg)

ell_radian, bee_radian, _, _ = \
ICRS_to_GAL(tmp_ra_radian,
            tmp_dec_radian,
            numpy.zeros_like(tmp_ra_radian),
            numpy.zeros_like(tmp_dec_radian) )


#################################################################
# select a subset of table by using "filter_good".
#
filter_good = (tmp_feh<-1) * (tmp_rv_err>0.) * (tmp_rv_err<20.)

data1 = data0[filter_good] #filter_good

col_ell_deg       = Column(name='l',        data = radian2deg(ell_radian[filter_good]) ) #filter_good
col_bee_deg       = Column(name='b',        data = radian2deg(bee_radian[filter_good]) ) #filter_good
col_HRV_corrected = Column(name='HRV_c',    data = data0['rv']  [filter_good] +  4.5   ) #filter_good
col_mean_mag      = Column(name='mean_mag', data = weighted_mean[filter_good]          ) #filter_good



table1 = Table(data1)

table1.add_columns([col_ell_deg,
                    col_bee_deg,
                    col_HRV_corrected,
                    col_mean_mag ])

table1['rv'].name = 'rv_original' # Fool-proof so that 'rv' cannot be mistakenly used.
table1['ra'].name = 'ra_LAMOST'   # Fool-proof to prevent confusion with 'ra' from Gaia DR2.
table1['dec'].name= 'dec_LAMOST'  # Fool-proof to prevent confusion with 'dec' from Gaia DR2.

# Sort table according to the value of 'rv_err'
# so that "unique(table1, keys='designation')" will only select the entry with lowest 'rv_error'.
table1.sort('rv_err')
table2 = unique(table1, keys='designation') # Select the entry with lowest 'rv_error'.

table2.write('LAMOSTDR5v3_feh-1_0rverr_RVcorrected_unique_meanmag.fits', overwrite=True)


```



---

```sql 

"# ADQL script run on TOPCAT"

SELECT 
mydb.*,
g.source_id, 
g.ra, 
g.dec, 
g.parallax, 
g.parallax_error,
g.pmra, 
g.pmra_error,
g.pmdec,
g.pmdec_error,
g.pmra_pmdec_corr,
g.phot_g_mean_mag, 
g.phot_bp_mean_mag, 
g.phot_rp_mean_mag,
g.radial_velocity,
g.radial_velocity_error,
distance(
  point('', mydb.ra_obs, mydb.dec_obs),
  point('', g.ra, g.dec)
) * 3600 as dist_arcsec
from TAP_UPLOAD.t55 as mydb 
join gaiadr2.gaia_source as g
  on 1 = contains(
    point('', mydb.ra_obs, mydb.dec_obs),
    circle('', g.ra, g.dec, 1. / 3600)
  ) 

```
===





## Some random thoughts.


- It is interesting to look at the metal-rich stellar substructure. For example, looking at retrograde, metal-rich stars. 

-