#! /usr/bin/env python
# -*- coding:utf8 -*-

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
from math import *

plt.style.use('classic')

# download data from ftp://ftp.marine.csiro.au/pub/legresy/gmsl_files/CSIRO_Alt.csv
data = []
with open('CSIRO_Alt.csv', 'r') as csv:
    for txt in csv.readlines():
        try:
            line = [i.strip() for i in txt.split(',')]
            if len(line) == 3:
                data.append(line)
        except Exception:
            pass

years0, h0 = [], []
years1, h1 = [], []
for d in data:
    try:
        y, h = float(d[0]), float(d[1])
        years0.append(y); h0.append(h)
    except Exception:
        pass
    try:
        y, h = float(d[0]), float(d[2])
        years1.append(y); h1.append(h)
    except Exception:
        pass

years0 = np.array(years0)

slope, intercept, r_value, p_value, std_err = stats.linregress(years0, h0)

plt.rc('font', size=15)
fig = plt.figure(figsize=(800 / 90.0, 540 / 90.0), dpi=72)
plt.plot(years0, h0, '-', color='#00eeee', linewidth=2, label='Monthly')
plt.plot(years1, h1, '-', color='#0000dd', linewidth=3, label='3-month running mean')
plt.plot(years0[[0,-1]], intercept+slope*np.array(years0[[0,-1]]), '-',
         color='red', linewidth=1.5, label='Trend = {:.1f} mm/year'.format(slope))
plt.grid(True)
plt.xlim(1992, ceil(max(years0)/1.)*1.)
plt.gca().ticklabel_format(useOffset=False)
plt.ylim(-50, 59)
plt.gca().ticklabel_format(useOffset=False)
plt.xlabel('Year')
plt.ylabel('Global Mean Sea Level (mm)')
plt.title('GMSL from TOPEX/Poseidon, Jason-1, Jason-2 and Jason-3 satellite altimeter data', size=14)
plt.legend(loc='upper left', prop={'size':14}, borderaxespad=1.2)
plt.text(0.47, 0.058, '''Seasonal signal removed
Inverse barometer correction applied
GIA correction applied''', size=14, transform=plt.gca().transAxes)
plt.tight_layout()
plt.savefig('Alt_gmsl_seas_rem.svg')