#!/usr/bin/env python

"""
A routine to calculate annual average values of mean sea level change using
linear regression.
"""

import numpy as np
import os
import glob
from sklearn.linear_model import LinearRegression

__author__ = "Ben Loveday"
__copyright__ = "Copyright 2023, EUMETSAT"
__credits__ = ["Ben Loveday", "Hayley Evers-King", "Remko Scharroo",
               "Stephen Killick"]
__license__ = "MIT"
__version__ = "1.0.1"
__maintainer__ = "Ben Loveday"
__email__ = "ops@eumetsat.int"
__status__ = "Production"

def make_averages(input_file, year_zero=2000, start_year=1993, end_year=2023):
    """
    A function to calculate annual average values of mean sea level change using
    linear regression.

    input_file (str) : path to the multi-mission time series
    year_zero (int)  : reference year, from data file header (default 2000)
    start_year (int) : analysis start year (default 1993)
    end_year (int)   : analysis end year (default 2023)

    return: None
    """

    # read the input file
    combined_record = np.genfromtxt(input_file, comments="#", dtype=float)

    # extract the times and global mean sea level
    times = np.array(combined_record[:, 0])
    MSL = np.array(combined_record[:, 1])

    # initialise lists
    years = []
    annual_rises = []

    # calculate annual mean sea level change values;
    for ii in np.arange(start_year - year_zero, end_year - year_zero):
        # find value in scope for this year
        kk = np.where((times > ii) & (times < ii+1))[0]
        # set the regression model
        model = LinearRegression().fit(times[kk].reshape((-1, 1)), MSL[kk])
        # extract the gradient and append to the list
        annual_rises.append(model.coef_)
        # append the year to the list
        years.append(year_zero + ii)

    # print the output
    for year, annual_rise in zip(years, annual_rises):
        print(f"{year}, {annual_rise[0]}")

if __name__ == '__main__':

    # set the input file
    input_file = "multimission_MSL.csv"

    # calculate the annual averages
    make_averages(input_file)