def main():
    #filenames = !ls *align/*reference-DEM.tif
    #run ~/src/demtools/error_analysis.py $filenames.s

    if len(sys.argv) < 1:
        sys.exit('No input files provided')

    fn_list = sys.argv[1:]
    n_samp = len(fn_list)

    error_dict_list = []
    for fn in fn_list:
        ed = parse_pc_align_log(fn)
        if 'Translation vector (North-East-Down, meters)' in ed.keys(): 
            error_dict_list.append(ed)  

    import matplotlib.dates as mdates
    #This is used for interactive display of x-value in plot window 
    date_str = '%Y-%m-%d %H:%M'
    date_fmt = mdates.DateFormatter(date_str)
    #ax.fmt_xdata = mdates.DateFormatter(date_fmt)
    months = mdates.MonthLocator() 
    months_int = mdates.MonthLocator(interval=6)  # every n months 
    years = mdates.YearLocator()   # every year
    yearsFmt = mdates.DateFormatter('%Y')
    #ax.xaxis.set_major_formatter(yearsFmt)
    #ax.xaxis.set_major_locator(months_int3)

    print
    print "n:", len(error_dict_list) 

    """
    #ECEF translations
    #key = 'Translation vector (ECEF meters)'
    key = 'Translation vector (Cartesian, meters)'
    #key = 'Translation vector (meters)'
    val = np.array([e[key] for e in error_dict_list])
    #make_plot3d(val[:,0], val[:,1], val[:,2], title=key)
    ce90 = geolib.CE90(val[:,0], val[:,1])
    le90 = geolib.LE90(val[:,2])
    print
    print key
    print "CE90:", ce90 
    print "LE90:", le90 
    print

    #Proj translation
    key = 'Translation vector (Proj meters)' 
    val = np.array([e[key] for e in error_dict_list])
    ce90 = geolib.CE90(val[:,0], val[:,1])
    le90 = geolib.LE90(val[:,2])
    print
    print key
    print "CE90:", ce90 
    print "LE90:", le90 
    print
    print 'Centroid (mean) of offsets (Proj meters): ', np.mean(val, axis=0) 
    print 'Centroid (median) of offsets (Proj meters): ', np.median(val, axis=0) 
    """

    #NOTE: changed default to N-E-D on 9/18/15
    #Can have significant differences for local proj vs. polar stereographic proj
    #Should regenerate all previous figures

    #Local translation on ellipsoid
    #This appears to be local stereographic projection on ellipsoid
    key = 'Translation vector (North-East-Down, meters)'
    val = np.array([e[key] for e in error_dict_list])

    #Reformat (n, e, +d) for (x, y, +z) coord sys
    val[:,[0,1]] = val[:,[1,0]]
    val[:,2] *= -1
    ce90 = geolib.CE90(val[:,0], val[:,1])
    le90 = geolib.LE90(val[:,2])
    print
    print key
    print "CE90:", ce90 
    print "LE90:", le90 
    print
    print 'Centroid (mean) of offsets (local ned meters): ', np.mean(val, axis=0) 
    print 'Centroid (median) of offsets (local ned meters): ', np.median(val, axis=0) 

    #Remove vertical bias
    remove_vertbias = False 
    if remove_vertbias:
        print "Removing vertical bias: %0.2f" % np.mean(val, axis=0)[2]
        val[:,2] -= np.mean(val, axis=0)[2]

    remove_outliers = False 
    #Flag outliers
    x_mag = val[:,0]
    y_mag = val[:,1]
    h_mag = np.sqrt(val[:,0]**2 + val[:,1]**2)
    v_mag = val[:,2]
    mag = np.sqrt(val[:,0]**2 + val[:,1]**2 + val[:,2]**2)
    abs_thresh = 10.0
    p = 98.0
    p_thresh = np.percentile(h_mag, p)
    #print "Outliers with horiz error >= %0.2f (%0.1f%%)" % (p_thresh, p)
    print "Outliers:" 
    #idx = (h_mag >= p_thresh).nonzero()[0]
    idx = (h_mag >= ce90).nonzero()[0]
    idx = np.unique(np.hstack([idx, ((np.abs(v_mag) >= le90).nonzero()[0])]))

    #Print all
    #idx = np.arange(h_mag.size)
    #idx_sort = np.argsort(mag[idx])
    #idx = idx[idx_sort]

    print 'name, m, h, v, x, y, z'
    for i in idx:
        print error_dict_list[i]['File'], mag[i], h_mag[i], v_mag[i], val[i,0:3]
        #Delete from list
        if remove_outliers:
            print "Removing from calculation"
            del error_dict_list[i]

    if remove_vertbias or remove_outliers:
        print
        print "Updated values"
        print key
        print "CE90:", geolib.CE90(val[:,0], val[:,1])
        print "LE90:", geolib.LE90(val[:,2])
        print
        print 'Centroid (mean) of offsets (local ned meters): ', np.mean(val, axis=0) 
        print 'Centroid (median) of offsets (local ned meters): ', np.median(val, axis=0) 

    #Extract dates
    date_vec = np.array([e['Date'] for e in error_dict_list])
    x = date_vec

    make_plot3d(val[:,0], val[:,1], val[:,2], title=key)
    #Note: there is a bug in pdf that displayes surface lines
    #fig_fn = 'icp_translation_vec_proj_meters.png'
    fig_fn = 'icp_translation_vec_local_meters.png'
    #plt.savefig(fig_fn, dpi=600, bbox_inches='tight')

    fig, ax = plt.subplots(1)
    key = 'Translation vector (lat,lon,z)'
    plt.title('ICP translation vector (lat,lon,z): Z component')
    val = np.array([e[key] for e in error_dict_list])
    y = val[:,2]
    make_plot(x,y,c='b',label=key, abs=False)
    fig.autofmt_xdate()
    ax.xaxis.set_minor_locator(months)
    #ax.xaxis.set_major_locator(months_int)
    #ax.xaxis.set_major_formatter(date_fmt)
    ax.fmt_xdata = date_fmt
    ax.set_ylabel('Z offset (m)')

    fig, ax = plt.subplots(1)
    key = 'Translation vector magnitude (meters)'
    plt.title('ICP Translation vector magnitude (meters)')
    y = np.array([e[key] for e in error_dict_list])
    make_plot(x,y,c='b',label=key, abs=True)
    fig.autofmt_xdate()
    ax.xaxis.set_minor_locator(months)
    #ax.xaxis.set_major_locator(months_int)
    #ax.xaxis.set_major_formatter(date_fmt)
    ax.fmt_xdata = date_fmt
    ax.set_ylabel('Offset (m)')

    fig, ax = plt.subplots(1)
    key = 'Number of errors'
    plt.title('Number of error samples')
    nerr = np.array([e[key] for e in error_dict_list])
    make_plot(x,nerr,c='b',label=key)
    fig.autofmt_xdate()
    ax.xaxis.set_minor_locator(months)
    #ax.xaxis.set_major_locator(months_int)
    #ax.xaxis.set_major_formatter(date_fmt)
    ax.fmt_xdata = date_fmt
    ax.set_ylabel('N samples')

    """
    fig, ax = plt.subplots(1)
    plt.title('ICP Standard Deviation')
    key = 'Input Std Error'
    in_std = np.array([e[key] for e in error_dict_list])
    make_plot(x,in_std,c='r',label=key)
    key = 'Output Std Error'
    out_std = np.array([e[key] for e in error_dict_list])
    make_plot(x,out_std,c='b',label=key)
    fig.autofmt_xdate()
    ax.xaxis.set_minor_locator(months)
    #ax.xaxis.set_major_locator(months_int)
    #ax.xaxis.set_major_formatter(date_fmt)
    ax.fmt_xdata = date_fmt
    plt.legend(scatterpoints=1)
    """

    fig, ax = plt.subplots(1)
    plt.title('ICP Mean Error')
    key = 'Input Mean Error'
    in_mean = np.array([e[key] for e in error_dict_list])
    #make_plot(x,in_mean,c='r',label=key,yerr=in_std)
    make_plot(x,in_mean,c='r',label=key, abs=True)
    key = 'Output Mean Error'
    out_mean = np.array([e[key] for e in error_dict_list])
    #make_plot(x,out_mean,c='b',label=key,yerr=out_std)
    make_plot(x,out_mean,c='b',label=key, abs=True)
    fig.autofmt_xdate()
    ax.xaxis.set_minor_locator(months)
    #ax.xaxis.set_major_locator(months_int)
    #ax.xaxis.set_major_formatter(date_fmt)
    ax.fmt_xdata = date_fmt
    ax.set_ylabel('Mean error (m)')
    plt.legend(scatterpoints=1, loc='upper left', prop={'size':8})

    fig, ax = plt.subplots(1)
    plt.title('ICP Median Error')
    key = 'Input 16th Percentile Error'
    in_16p = np.array([e[key] for e in error_dict_list])
    key = 'Input 84th Percentile Error'
    in_84p = np.array([e[key] for e in error_dict_list])
    key = 'Input Median Error'
    in_med = np.array([e[key] for e in error_dict_list])
    make_plot(x,in_med,c='r',label=key,yerr=[in_med - in_16p, in_84p - in_med], abs=True)
    key = 'Output 16th Percentile Error'
    out_16p = np.array([e[key] for e in error_dict_list])
    key = 'Output 84th Percentile Error'
    out_84p = np.array([e[key] for e in error_dict_list])
    key = 'Output Median Error'
    out_med = np.array([e[key] for e in error_dict_list])
    make_plot(x,out_med,c='b',label=key,yerr=[out_med - out_16p, out_84p - out_med], abs=True)
    fig.autofmt_xdate()
    ax.fmt_xdata = mdates.DateFormatter(date_fmt)
    ax.xaxis.set_minor_locator(months)
    #ax.xaxis.set_major_locator(months_int)
    #ax.xaxis.set_major_formatter(date_fmt)
    ax.fmt_xdata = date_fmt
    ax.set_ylabel('Median error (m)')
    plt.legend(scatterpoints=1, loc='upper left', prop={'size':8})
    fig_fn = 'icp_median_error.pdf'
    plt.savefig(fig_fn, dpi=600, bbox_inches='tight')

    fig, ax = plt.subplots(1)
    plt.title('Sampled Median Error')
    key = 'Input Sampled 16th Percentile Error'
    in_16p = np.ma.fix_invalid([e[key] for e in error_dict_list])
    if in_16p.count() > 0:
        key = 'Input Sampled 84th Percentile Error'
        in_84p = np.ma.fix_invalid([e[key] for e in error_dict_list])
        key = 'Input Sampled Median Error'
        in_med = np.ma.fix_invalid([e[key] for e in error_dict_list])
        in_spread = in_84p - in_16p
        make_plot(x,in_med,c='r',label=key,yerr=[in_med - in_16p, in_84p - in_med], abs=True)

        key = 'Output Sampled 16th Percentile Error'
        out_16p = np.ma.fix_invalid([e[key] for e in error_dict_list])
        key = 'Output Sampled 84th Percentile Error'
        out_84p = np.ma.fix_invalid([e[key] for e in error_dict_list])
        key = 'Output Sampled Median Error'
        out_med = np.ma.fix_invalid([e[key] for e in error_dict_list])
        out_spread = out_84p - out_16p

        p = 90.0
        out_med_thresh = np.percentile(out_med, p)
        out_spread_thresh = np.percentile(out_spread, p)
        #print "Outliers with horiz error >= %0.2f (%0.1f%%)" % (p_thresh, p)
        print
        print "Sampled Error Outliers:" 
        #idx = (h_mag >= p_thresh).nonzero()[0]
        idx = (out_med >= out_med_thresh).nonzero()[0]
        idx = np.unique(np.hstack([idx, ((out_spread >= out_spread_thresh).nonzero()[0])]))
        #Print all
        idx = np.arange(out_med.size)
        idx_sort = np.argsort(out_med[idx])
        idx = idx[idx_sort]
        print 'name, samp_mederrr, samp_errspread, nerr'
        for i in idx:
            print error_dict_list[i]['File'], out_med[i], out_spread[i], nerr[i]
            #Delete from list
            if remove_outliers:
                print "Removing from calculation"
                del error_dict_list[i]
        print
        print 'Input sampled median error (spread/2): %0.2f (%0.2f)' % (np.median(in_med), np.median(in_spread)/2.)
        print 'Output sampled median error (spread/2): %0.2f (%0.2f)' % (np.median(out_med), np.median(out_spread)/2.)
        print

        make_plot(x,out_med,c='b',label=key,yerr=[out_med - out_16p, out_84p - out_med], abs=True)
        fig.autofmt_xdate()
        ax.set_ylabel('Median error (m)')
        ax.fmt_xdata = mdates.DateFormatter(date_fmt)
        ax.xaxis.set_minor_locator(months)
        #ax.xaxis.set_major_locator(months_int)
        #ax.xaxis.set_major_formatter(date_fmt)
        ax.fmt_xdata = date_fmt
        ax.set_ylabel('Median error (m)')
        plt.legend(scatterpoints=1, loc='upper left', prop={'size':8})
        ax.set_ylim(-15,15)
        fig_fn = 'sampled_median_error.pdf'
        #fig_fn = 'sampled_median_error_2014-2016.pdf'
        #from datetime import datetime
        #ax.set_xlim(datetime(2014,1,1),datetime(2016,7,1))
        plt.savefig(fig_fn, dpi=600, bbox_inches='tight')
#ECEF translations
    #key = 'Translation vector (ECEF meters)'
    key = 'Translation vector (Cartesian, meters)'
    #key = 'Translation vector (meters)'
    val = np.array([e[key] for e in error_dict_list])
    #make_plot3d(val[:,0], val[:,1], val[:,2], title=key)
    ce90 = geolib.CE90(val[:,0], val[:,1])
    le90 = geolib.LE90(val[:,2])
    print
    print key
    print "CE90:", ce90 
    print "LE90:", le90 
    print

    #Proj translation
    key = 'Translation vector (Proj meters)' 
    val = np.array([e[key] for e in error_dict_list])
    ce90 = geolib.CE90(val[:,0], val[:,1])
    le90 = geolib.LE90(val[:,2])
    print
    print key
    print "CE90:", ce90 
    print "LE90:", le90 
    print
    print 'Centroid (mean) of offsets (Proj meters): ', np.mean(val, axis=0) 
    print 'Centroid (median) of offsets (Proj meters): ', np.median(val, axis=0)