### St John Wind Speed Processing Code
# Alexander Rey, June 30, 2020

from erddapy import ERDDAP
import pandas as pd
import datetime
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import scipy as sc
import statsmodels.api as sm

# Import Offshore Buoy Data
e = ERDDAP(
    server='https://www.smartatlantic.ca/erddap',
    protocol='tabledap',
    response='csv',
)

# Setup download for the buoy
e.dataset_id = 'SMA_saint_john'
e.constraints = {
    'time>=': '2016-07-10T00:00:00Z',
    'time<=': '2020-06-01T00:00:00Z',
}
e.variables = [
    'time',
    'wind_spd_avg',
    'wind_dir_avg',
]

buoydat = e.to_pandas(
    index_col='time (UTC)',
    parse_dates=True,
).dropna()

# Setup download for the warf
e.dataset_id = 'SMA_saint_john_wharf'
e.constraints = {
    'time>=': '2016-07-10T00:00:00Z',
    'time<=': '2020-06-01T00:00:00Z',
}
e.variables = [
    'time',
    'wind_spd_avg',
    'wind_dir_avg',
]

warfdat = e.to_pandas(
    index_col='time (UTC)',
    parse_dates=True,
).dropna()

# Setup download for the cruise terminal
e.dataset_id = 'SMA_saint_john_cruise_terminal'
e.constraints = {
    'time>=': '2016-12-23T00:00:00Z',
    'time<=': '2020-06-01T00:00:00Z',
}

e.variables = [
    'time',
    'wind_spd_avg',
    'wind_dir_avg',
]

cruisedat = e.to_pandas(
    index_col='time (UTC)',
    parse_dates=True,
).dropna()

### Clean and process data
t_index = pd.date_range(start='2016-07-10T00:00:00Z', end='2020-06-01T00:00:00Z', freq='30Min')

# Set to 30 minute time series
buoy_resample = buoydat.resample('30Min').mean().reindex(t_index)
buoy_resample['wind_spd_avg (m s-1)'] = buoy_resample['wind_spd_avg (m s-1)'].interpolate('linear')
buoy_resample['wind_dir_avg (degree)'] = buoy_resample['wind_dir_avg (degree)'].interpolate('nearest')

warf_resample = warfdat.resample('30Min').mean().reindex(t_index)
warf_resample['wind_spd_avg (m s-1)'] = warf_resample['wind_spd_avg (m s-1)'].interpolate('linear')
warf_resample['wind_dir_avg (degree)'] = warf_resample['wind_dir_avg (degree)'].interpolate('nearest')

cruise_resample = cruisedat.resample('30Min').mean().reindex(t_index)
cruise_resample['wind_spd_avg (m s-1)'] = cruise_resample['wind_spd_avg (m s-1)'].interpolate('linear')
cruise_resample['wind_dir_avg (degree)'] = cruise_resample['wind_dir_avg (degree)'].interpolate('nearest')




# %% Make QQ plot
plt.figure()
# plt.scatter(buoy_resample['wind_spd_avg (m s-1)'], warf_resample['wind_spd_avg (m s-1)'],
#            1,buoy_resample['wind_dir_avg (degree)'],marker='.')

# sm.qqplot_2samples(warf_resample['wind_spd_avg (m s-1)'].loc[buoy_resample['wind_dir_avg (degree)'].isin([150, 220])],
#                    buoy_resample['wind_spd_avg (m s-1)'].loc[buoy_resample['wind_dir_avg (degree)'].isin([150, 220])],
#                    line='r',ylabel='Potash Terminal Quantiles',xlabel='Offshore Quantiles')
# plt.title('St John Wind (Potash Terminal) from 160-220 deg July 2016-June 2020')


sm.qqplot_2samples(cruise_resample['wind_spd_avg (m s-1)'].loc[buoy_resample['wind_dir_avg (degree)'].isin([160, 220])],
                   buoy_resample['wind_spd_avg (m s-1)'].loc[buoy_resample['wind_dir_avg (degree)'].isin([160, 220])],
                   line='r', ylabel='Cruise Terminal Quantiles', xlabel='Offshore Quantiles')
plt.title('St John Wind (Cruise Terminal) from 160-220 deg Dec 2016-June 2020')

# plt.scatter(buoy_resample['wind_spd_avg (m s-1)'].loc[buoy_resample['wind_dir_avg (degree)'].isin([90, 180])],
#             warf_resample['wind_spd_avg (m s-1)'].loc[buoy_resample['wind_dir_avg (degree)'].isin([90, 180])],
#             10, marker='.')

plt.plot([0, 20], [0, 20], 'k-')

# plt.xlabel('Offshore Wind Speed (m/s)')
# plt.ylabel('Potash Terminal Wind Speed (m/s)')
# plt.title('St John Wind Speed Comparision July 2016-June 2020')
# plt.title('St John Wind Speed Comparision July 2016-June 2020, 90-180 deg')

plt.xlim(0, 20)
plt.ylim(0, 20)

# cmap = plt.get_cmap('hsv')
# norm = plt.Normalize(0, 360)
# color = cmap(norm(200.))
#
# cbar = plt.colorbar()
# cbar.ax.get_yaxis().labelpad = 15
# cbar.ax.set_ylabel('Offshore Wind Direction', rotation=270)

plt.savefig('St_John_Cruise_160_220.png')
# plt.savefig('St_John_all.png')
plt.show()