data_filterer.py

"""
Script to be used for filtering of corrected flux data files generated by datasorter.py.
- Loads corrected flux and time from file.
- For each data-point, calculates a median flux value from the surrounding 9 points (including itself, 4 to each side).
- Plots variance of data from median value.
- Optional: Allows you to graphically select a specific region to exclude from the data-set.
- Optional: Can automatically delete data that is surrounded by bad data points. Not recommended, not fully tested.
- Calculates standard deviation of variance from median.
- After this, will automatically delete data points that has a variance larger than n standard deviations (n chosen by
  text prompt).
- Saves filtered data to file.
"""

import numpy as np
import matplotlib.pyplot as plt
import ps_f
from scipy.ndimage import median_filter

# Plot parameters (font size and option to skip plotting)
plt.rcParams.update({'font.size': 22})  # Changes font size on all plots
skip_plot = False

# # Data reading
inpt1 = input('What is the filename (without extension)?')  # Input for data reader
data = ps_f.reader(inpt1)

time = data[0]
flux = data[1]
flux = flux - np.mean(flux)

start = time[0]
daytype = 'JD'
time = time - start  # To set time[0]=0 as reference
time_days = np.copy(time)
time = time * 86400  # change to time in seconds


# # Median points
mflux = median_filter(flux, 9, mode='reflect')

# # Median plot
if skip_plot is not True:
    plt.figure()
    plt.plot(time_days, flux, 'b.')
    plt.plot(time_days, flux, 'b--', linewidth=0.5)
    plt.xlabel('Days after ' + str(start) + daytype)
    plt.ylabel('Corrected flux ppm')
    plt.plot(time_days, mflux, 'g.')
    plt.legend(['Data', 'Median of data from surrounding points'])
    plt.show(block=False)


# #  Distribution of points

variance = flux - mflux
h = np.sort(variance)  # Variance sorted from low to high value


# # Plot variance
if skip_plot is not True:
    plt.figure()
    m = len(variance)
    plt.plot(time_days, variance, 'r.')
    plt.xlabel('Days after ' + str(start) + daytype)
    plt.ylabel('Variance')
    plt.show(block=False)

    plt.figure()
    plt.plot(sorted(variance), 'r.')
    plt.ylabel('sorted variance')
    plt.show(block=True)


# # Deletion of data points by manual selection of interval
print('Do you want to remove points in a region manually? (y/n)')
inpt4 = input()
if inpt4 == 'y':
    # # Graphical input selection of region to delete
    plt.figure()
    plt.plot(variance, 'b.')
    plt.xlabel('List element')
    plt.ylabel('Variance')
    plt.title('Select region which will be deleted (only 1 region, two points)')
    coords = plt.ginput(n=2, timeout=0, show_clicks=True, mouse_add=1, mouse_stop=3, mouse_pop=2)  # graphical input

    # # Conversion of from input coordinates to x-axis integer coordinates, x-axis is list index
    [coord1, coord2] = [coords[0], coords[1]]
    coord1 = int(coord1[0])
    coord2 = int(coord2[0])

    del_select = list(range(coord1, coord2))  # List with deletion indices to be used for plotting

    # # Plot selected region
    if skip_plot is not True:
        plt.figure()
        plt.plot(variance, 'r.')
        plt.plot(del_select, variance[coord1:coord2], 'g.')
        plt.show(block=False)

    # # Delete selected region from data set
    variance = np.delete(variance, del_select)
    time = np.delete(time, del_select)
    time_days = np.delete(time_days, del_select)
    flux = np.delete(flux, del_select)
    mflux = np.delete(mflux, del_select)

    if skip_plot is not True:
        # # Plotting of variance after deletion of data
        plt.figure()
        plt.plot(time_days, variance, 'r.')
        plt.xlabel('Days after ' + str(start) + daytype)
        plt.ylabel('Variance')
        plt.show(block=False)

        # # Plotting of sorted variance after deletion of data
        plt.figure()
        plt.plot(sorted(variance), 'r.')
        plt.ylabel('sorted variance')
        plt.show()


# # Input to decide what the upper limit for deleting a point is (indicated in stand dev, usually 3 or 4 is used)
inpt2 = float(input('How many standard deviations away do we short out data?'))  # Must be an integer


# #  Deletion of points further than inpt2 std away in variance
stdev = np.std(variance)  # Standard deviation of variance

del_select = np.where(np.abs(variance) > inpt2 * stdev)[0]

variance = np.delete(variance, del_select)
time = np.delete(time, del_select)
time_days = np.delete(time_days, del_select)
flux = np.delete(flux, del_select)
mflux = np.delete(mflux, del_select)


if skip_plot is not True:
    # # Plot new sorted variance
    plt.figure()
    plt.plot(sorted(variance), 'r.')
    plt.ylabel('sorted variance after discarding points')
    plt.show(block=False)

    # # Plot new variance vs. time
    plt.figure()
    plt.plot(time_days, variance, 'r.')
    plt.xlabel('Days after ' + str(start) + daytype)
    plt.ylabel('Variance')
    plt.show(block=False)

    # # Plot new data set (flux)
    plt.figure()
    plt.plot(time_days, flux, 'b.')
    plt.plot(time_days, mflux, 'g.')
    plt.xlabel('Days after ' + str(start) + daytype)
    plt.ylabel('Corrected flux (ppm)')
    plt.show(block=True)

# # Writer call to save data
inpt3 = input('What is the filename (without extension)?')
plt.close('all')
ps_f.writer(inpt3, time, flux, mflux, variance)