from obspy import read as rdseed
import numpy as np
from scipy.stats import skew, kurtosis
from scipy import fftpack
# import pandas as pd

def calc_features(mseed_file, start_time, end_time, jj):

    # Reading waveforms

    continuous_stream = rdseed(mseed_file, type='mseed')
    continuous_stream.filter('bandpass', freqmin=1.0, freqmax=15.0, corners=2, zerophase=True)
    str_fs = float(continuous_stream.traces[0].stats.sampling_rate)

    # Trimming waveforms
    stream_trimmed = continuous_stream.trim(start_time, end_time)
    data_trimmed = stream_trimmed.traces[0].data
    data_final = data_trimmed.copy()

    # if jj == 466:
    #     stream_trimmed.plot()

    # Time Domain Features

    # Feature 1
    dur_array_t = (len(data_final) - 1) / float(str_fs)

    # Feature 2
    avg_array_t = np.mean(data_final)

    # Feature 3
    std_array_t = np.std(data_final)

    # Feature 4
    skw_array_t = skew(data_final)

    # Feature 5
    kurt_array_t = kurtosis(data_final)

    # Feature 6
    idx_max_t = np.where(data_final == np.max(data_final))[0]
    ratio_ascdesc = (float(idx_max_t) / str_fs) / (dur_array_t - (float(idx_max_t) / str_fs))


    time_feat = np.array([dur_array_t, avg_array_t, std_array_t, skw_array_t, kurt_array_t,
                          ratio_ascdesc])
    time_features = time_feat[np.newaxis, :]



    # Frequency Domain
    y_dft = fftpack.fft(data_final)
    x_freq = fftpack.fftfreq(data_final.size, d=(1 / str_fs))
    power_amp = (np.abs(y_dft)) / len(data_final)


    # Feature 7
    id_maxf = np.where(power_amp == np.max(power_amp))[0]
    dom_freq = x_freq[id_maxf][0]

    # Feature 8
    max_spect = np.max(power_amp)
    median_norm = np.median(power_amp/max_spect)

    freq_feat = np.array([dom_freq, median_norm])
    freq_features = freq_feat[np.newaxis, :]

    tf_features = np.concatenate([time_features, freq_features], axis=1)

    return tf_features