Mm rework lingfeng

.gitignore

@@ -1,5 +1,6 @@
 *.png
 *random*.fits
+*random*.lis
 *.DS_Store
 *.pyc
 .ipynb_checkpoints/

flystar/analysis.py

@@ -1,17 +1,11 @@
 import numpy as np
 import pylab as plt
-from flystar import starlists
-from flystar import startables
-from flystar import align
-from flystar import match
-from flystar import transforms
+from . import starlists, match
 from astropy import table
 from astropy.table import Table, Column
 from astropy.coordinates import SkyCoord
 from astropy import units as u
 from astropy.wcs import WCS
-from astroquery.gaia import Gaia
-from astroquery.mast import Observations, Catalogs
 import pdb, copy
 import math
 from scipy.stats import f
@@ -35,25 +29,27 @@ def query_gaia(ra, dec, search_radius=30.0, table_name='gaiadr3'):
         Dec. in degrees in the format such as '-29:00:28.0'
 
     search_radius : float
-        The search radius in arcseconds. 
+        The search radius in arcseconds.
 
     Optional Input
     --------------
     table_name : string
         Options are 'gaiadr2' or 'gaiaedr3'
     """
+    from astroquery.gaia import Gaia
     target_coords = SkyCoord(ra, dec, unit=(u.hourangle, u.deg), frame='icrs')
     ra = target_coords.ra.degree
     dec = target_coords.dec.degree
 
     search_radius *= u.arcsec
 
     Gaia.ROW_LIMIT = 50000
-    gaia_job = Gaia.cone_search_async(target_coords, search_radius, table_name = table_name + '.gaia_source')
+    gaia_job = Gaia.cone_search_async(target_coords, radius=search_radius, table_name=table_name + '.gaia_source')
     gaia = gaia_job.get_results()
 
     #Change new 'SOURCE_ID' column header back to lowercase 'source_id' so all subsequent functions still work:
-    gaia['SOURCE_ID'].name = 'source_id'
+    if 'SOURCE_ID' in gaia.colnames:
+        gaia.rename_column('SOURCE_ID', 'source_id')
 
     return gaia
 
@@ -107,13 +103,13 @@ def check_gaia_parallaxes(ra,dec,search_radius=10.0,table_name='gaiadr3',target=
     plt.yscale('log')
     plt.tight_layout()
     plt.savefig('gaiaplx'+file_ext+'.png')
-    
+
 
 def prepare_gaia_for_flystar(gaia, ra, dec, targets_dict=None, match_dr_max=0.2, pi_err_limit=0.4, default_motion_model='Linear'):
     """
     Take a Gaia table (from astroquery) and produce a new table with a tangential projection
-    and shift such that the origin is centered on the target of interest. 
-    Convert everything into arcseconds and name columns such that they are 
+    and shift such that the origin is centered on the target of interest.
+    Convert everything into arcseconds and name columns such that they are
     ready for FlyStar input.
 
     Inputs
@@ -130,7 +126,7 @@ def prepare_gaia_for_flystar(gaia, ra, dec, targets_dict=None, match_dr_max=0.2,
     target_coords = SkyCoord(ra, dec, unit=(u.hourangle, u.deg), frame='icrs')
     ra = target_coords.ra.degree     # in decimal degrees
     dec = target_coords.dec.degree   # in decimal degrees
-    
+
     cos_dec = np.cos(np.radians(dec))
     x = (gaia['ra'] - ra) * cos_dec * 3600.0   # arcsec
     y = (gaia['dec'] - dec) * 3600.0           # arcsec
@@ -149,7 +145,7 @@ def prepare_gaia_for_flystar(gaia, ra, dec, targets_dict=None, match_dr_max=0.2,
     gaia_new['vy'] = gaia['pmdec'].data / 1e3
     gaia_new['vx_err'] = gaia['pmra_error'].data / 1e3
     gaia_new['vy_err'] = gaia['pmdec_error'].data / 1e3
-    
+
     gaia_new['t0'] = gaia['ref_epoch'].data
     gaia_new['source_id'] = gaia['source_id'].data.astype('S19')
 
@@ -159,7 +155,7 @@ def prepare_gaia_for_flystar(gaia, ra, dec, targets_dict=None, match_dr_max=0.2,
     gaia_new['vy'][idx] = 0.0
     gaia_new['vx_err'][idx] = 0.0
     gaia_new['vy_err'][idx] = 0.0
-    
+
     gaia_new['m'] = gaia['phot_g_mean_mag']
     gaia_new['me'] = 1.09/gaia['phot_g_mean_flux_over_error']
     gaia_new['pi'] = gaia['parallax'].data*1e-3
@@ -171,7 +167,7 @@ def prepare_gaia_for_flystar(gaia, ra, dec, targets_dict=None, match_dr_max=0.2,
     gaia_new['vx_err'][idx] = 0.0
     gaia_new['vy'][idx] = 0.0
     gaia_new['vy_err'][idx] = 0.0
-    
+
     # Cut out stars with high plx error and set motion models
     idx = np.where((gaia_new['pi_err']>(pi_err_limit/1e3)) | (gaia['parallax'].mask == True))[0]
     gaia_new['pi'][idx] = 0.0
@@ -190,9 +186,13 @@ def prepare_gaia_for_flystar(gaia, ra, dec, targets_dict=None, match_dr_max=0.2,
         gaia_new['motion_model_input'] = 'Fixed'
         gaia_new['motion_model_used'] = 'Fixed'
         gaia_new['n_params'] = 1
+    elif default_motion_model=='Empty':
+        gaia_new['motion_model_input'] = 'Empty'
+        gaia_new['motion_model_used'] = 'Empty'
+        gaia_new['n_params'] = 0
     else:
         print("Invalid motion model",default_motion_model,"- none assigned")
-        
+
     #macy additions to try to fix wild magnitude values
     #gaia_new['ruwe'] = gaia['ruwe']
     #try:
@@ -228,9 +228,9 @@ def prepare_gaia_for_flystar(gaia, ra, dec, targets_dict=None, match_dr_max=0.2,
                 print('Found match for: ', targ_names[idx], ' - ',gaia_new['source_id'][i_gaia])
 
     return gaia_new
-    
+
 def run_flystar():
-    
+
     test_file = '/u/jlu/work/microlens/OB150211/a_2018_10_19/a_ob150211_2018_10_19/lis/stars_matched2.fits'
 
     t = Table.read(test_file)
@@ -262,39 +262,39 @@ def run_flystar():
     ym_t = y0 + vy * (t - t0)
 
     # Model distorted positions
-    
-    
+
+
     return
 
 
 def project_gaia(gaia, epoch, ra, dec):
     """
     Take the Gaia measurements, forward them in time, and then convert them into a tangential projection.
-    
+
     Inputs
     ----------
     epoch : float (year)
         The decimal year to project the measurement to. Note that we use 365.25 days per year.
-        
+
     ra : float (deg)
         The right ascension (J2000) in decimal degrees of the center of the field.
-        
+
     dec : float (deg)
         The declination (J2000) in decimal degrees of the center of the field.
-        
+
     """
     t0 = gaia['ref_epoch']
     x0 = (gaia['ra'] - ra) * np.cos(np.radians(dec)) * 3600.0   # Arcsec
     y0 = (gaia['dec'] - dec) * 3600.0
     x0e = gaia['ra_error']  / 1.0e3       # arcsec, already in alpha* (multiplied by cos(delta))
     y0e = gaia['dec_error'] / 1.0e3       # arcsec
-    
-    
+
+
     vx = gaia['pmra'] / 1.0e3            # arcsec / yr
-    vy = gaia['pmdec'] / 1.0e3       
+    vy = gaia['pmdec'] / 1.0e3
     vxe = gaia['pmra_error'] / 1.0e3     # arcsec / yr
     vye = gaia['pmdec_error'] / 1.0e3
-    
+
     # Modify any vx/vy, etc. that are zero and make a regular (unmasked) numpy array.
     vx[vx.mask] = 0.0
     vy[vy.mask] = 0.0
@@ -304,29 +304,29 @@ def project_gaia(gaia, epoch, ra, dec):
     vy = np.array(vy)
     vxe = np.array(vxe)
     vye = np.array(vye)
-    
+
     dt = epoch - t0
     x_now = (x0 + (vx * dt)) * -1.0  # Switch to a left-handed coordinate system, like detector pixels.
     y_now = (y0 + (vy * dt))
     xe_now = np.hypot(x0e, vxe*dt)
     ye_now = np.hypot(y0e, vye*dt)
-    
+
     # Format as a starlist
-    gaia_lis = starlists.StarList(name=gaia['source_id'], 
+    gaia_lis = starlists.StarList(name=gaia['source_id'],
                                   x=x_now, y=y_now, m=gaia['phot_g_mean_mag'],
                                   xe=xe_now, ye=ye_now, me=1.0/gaia['phot_g_mean_flux_over_error'])
-    
+
     # Duplicate columns to 'x_avg', etc. Needed for initial guessing.
     gaia_lis['x_avg'] = gaia_lis['x']
     gaia_lis['y_avg'] = gaia_lis['y']
-    gaia_lis['m_avg'] = gaia_lis['m']    
-    
+    gaia_lis['m_avg'] = gaia_lis['m']
+
     return gaia_lis
 
 
 def rename_after_flystar(star_tab, label_dat_file, new_copy=True, dr_tol=0.05, dm_tol=0.3, verbose=False):
     """
-    Take a StarTable output from FlyStar MosaicToRef that has been 
+    Take a StarTable output from FlyStar MosaicToRef that has been
     aligned into R.A. and Dec. (usually by way of Gaia). Align
     the output to a label.dat file for this source and rename
     everything.
@@ -354,28 +354,28 @@ def rename_after_flystar(star_tab, label_dat_file, new_copy=True, dr_tol=0.05, d
                              x_lab[ndx_lab[ii]], star_tab['x0'][ndx_star[ii]],
                              y_lab[ndx_lab[ii]], star_tab['y0'][ndx_star[ii]],
                              m_lab[ndx_lab[ii]], star_tab['m0'][ndx_star[ii]]))
-        
+
 
         print('Temporary shift transformations: ')
         print('    dm = {0:8.4f} +/- {1:8.4f}'.format(dm.mean(), dm.std()))
         print('    dx = {0:8.4f} +/- {1:8.4f}'.format(dx.mean(), dx.std()))
         print('    dy = {0:8.4f} +/- {1:8.4f}'.format(dy.mean(), dy.std()))
-    
+
     m_lab = label_tab['m'] + dm.mean()
     x_lab += dx.mean()
     y_lab += dy.mean()
-    
+
     # Now that we are in a common coordinate and magnitude
     # system, lets match the whole lists by coordinates.
-    idx_lab, idx_star, dr, dm = match.match(x_lab, y_lab, m_lab, 
+    idx_lab, idx_star, dr, dm = match.match(x_lab, y_lab, m_lab,
                                             star_tab['x0'], star_tab['y0'], star_tab['m0'],
                                             dr_tol=dr_tol, dm_tol=dm_tol, verbose=verbose)
     #print('idx_lab:')
     #for iii in range(len(idx_lab)):
     #    print(label_tab["name"][idx_lab[iii]], star_tab["name"][idx_star[iii]])
 
     print('Renaming {0:d} out of {1:d} stars'.format(len(idx_lab), len(star_tab)))
-    
+
     # Make a copy of the table, UNLESS, the user specifies.
     if new_copy:
         new_tab = copy.deepcopy(star_tab)
@@ -385,9 +385,9 @@ def rename_after_flystar(star_tab, label_dat_file, new_copy=True, dr_tol=0.05, d
     # copy over the original names... don't overwrite (this could mean data loss)
     if 'name_orig' not in new_tab.colnames:
         new_tab.add_column(Column(star_tab['name'].data, name='name_orig'))
-        
+
     new_tab['name'][idx_star] = label_tab[idx_lab]['name']
-    
+
     return new_tab
 
 def pick_good_ref_stars(star_tab, r_cut=None, m_cut=None, p_err_cut=None, pm_err_cut=None, name_cut=None, reset=True):
@@ -432,9 +432,9 @@ def pick_good_ref_stars(star_tab, r_cut=None, m_cut=None, p_err_cut=None, pm_err
 
 def startable_subset(tab, idx, mag_trans=True, mag_trans_orig=False):
     """
-    Input is MosaicToRef table from alignment of multiple filters, 
+    Input is MosaicToRef table from alignment of multiple filters,
     such that the astrometry is combined but the photometry is not.
-    This function is used to separate out a selected filter from the 
+    This function is used to separate out a selected filter from the
     combined astrometry + uncombined photometry table.
     """
     # Multiples: ['x', 'y', 'm', 'name_in_list', 'xe', 'ye', 'me', 't',
@@ -466,7 +466,7 @@ def startable_subset(tab, idx, mag_trans=True, mag_trans_orig=False):
             # Update the original table.
             if mag_trans_orig:
                 tab['m'][:,idx[ii]] += mag_offset
-    
+
     return new_tab
 
 
@@ -476,13 +476,13 @@ def startable_subset(tab, idx, mag_trans=True, mag_trans_orig=False):
 
 def calc_chi2(ref_mat, starlist_mat, transform, errs='both'):
     """
-    calculate the chi2 and reduced chi2 of the position 
+    calculate the chi2 and reduced chi2 of the position
     between two matched starlists.
     Input:
     ref_mat: astropy table
         Reference starlist only containing matched stars that were used in the
         transformation. Standard column headers are assumed.
-        
+
     starlist_mat: astropy table
         Transformed starlist only containing the matched stars used in
         the transformation. Standard column headers are assumed.
@@ -522,19 +522,19 @@ def calc_chi2(ref_mat, starlist_mat, transform, errs='both'):
     elif errs == 'starlist':
         xerr = starlist_mat['xe']
         yerr = starlist_mat['ye']
-          
+
 
     # For both X and Y, calculate chi-square. Combine arrays to get combined
     # chi-square
     chi_sq_x = diff_x**2. / xerr**2.
     chi_sq_y = diff_y**2. / yerr**2.
 
     chi_sq = np.append(chi_sq_x, chi_sq_y)
-    
+
     # Calculate degrees of freedom in transformation
     num_mod_params = calc_nparam(transform)
     deg_freedom = len(chi_sq) - num_mod_params
-    
+
     # Calculate reduced chi-square
     chi_sq = np.sum(chi_sq)
     chi_sq_red = chi_sq / deg_freedom
@@ -551,7 +551,7 @@ def calc_nparam(transformation):
         nparam = 4
     elif transformation.__class__.__name__ == 'PolyTransform':
         order = transformation.order
-        nparam = (order+1) * (order+2) 
+        nparam = (order+1) * (order+2)
     return nparam
 
 def calc_F(red_chi2_1, red_chi2_2, v1, v2):
@@ -572,24 +572,24 @@ def calc_F(red_chi2_1, red_chi2_2, v1, v2):
     for 1st order polynomial fitting:
         x' = a0 + a1*x + a2*y
         y' = b0 + b1*x + b2*y
-        v1 = 2*N1 - 2*3 (2*: because x and y direction) 
+        v1 = 2*N1 - 2*3 (2*: because x and y direction)
         red_chi2_1 = chi2/v1
     for 2nd order polynomial fitting:
         x' = a0 + a1*x + a2*y + a3*x**2 + a4*y**2 + a5*x*y
         y' = b0 + b1*x + b2*y + b3*x**2 + b4*y**2 + b5*x*y
-        v1 = 2*N1 - 2*6 
+        v1 = 2*N1 - 2*6
         red_chi2_2 = chi2/v2
     calc_F(red_chi2_1, red_chi2_2, v1, v2)
-    
+
     ***Note***
-    * make sure the first model is the simple model 
+    * make sure the first model is the simple model
       and the second model is the more complicated model
-    * the return value represents the probability that 
+    * the return value represents the probability that
       the first model is better than the second model, in other words,
       the small P means the more colicated model is needed.
       the large P means the simple model is good enough.
-    * normally, the P value will increase from model1->model2, to 
-      model2->model3, to model3->model4. The user can decide a 
+    * normally, the P value will increase from model1->model2, to
+      model2->model3, to model3->model4. The user can decide a
       critical value (eg, 0.7) to find the proper model.
     """
 

flystar/archive_io.py

@@ -1,17 +1,17 @@
 import pickle
 
-# Need to add these functions to a utility .py file rather than storing them in general structure. 
+# Need to add these functions to a utility .py file rather than storing them in general structure.
 def open_archive(file_name):
     """
-    Helper function to open archived files. 
+    Helper function to open archived files.
     """
     with open(file_name, 'rb') as file_archive:
         file_dict = pickle.load(file_archive)
     return file_dict
 
 def save_archive(file_name, save_data):
     """
-    Helper function to archive a file. 
+    Helper function to archive a file.
     """
     with open(file_name, 'wb') as outfile:
         pickle.dump(save_data, outfile, protocol=pickle.HIGHEST_PROTOCOL)