# Notebook setup
import xarray as xr
!pip install netcdf4
import matplotlib.pyplot as plt

# This makes the plots prettier
import seaborn as sns
sns.set()
sns.set(rc={'figure.figsize':(8,6)})

# This removes an annoying warning
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

# The OPENDAP URL for the station we want
url = 'https://dods.ndbc.noaa.gov/thredds/dodsC/data/stdmet/44025/44025.ncml'

# Open the dataset using xarray, look how simple!
ds = xr.open_dataset(url)

# Quick list of available variaibles
ds.data_vars

# Limit to last decade
ds = ds.sel(time=slice('2010-01-01', '2020-01-01'))

# Quickplot of SST
ds.sea_surface_temperature.plot()
plt.title('NDBC Buoy 44025');

# Calculate and Plot Daily Average
daily_sst = ds.sea_surface_temperature.load().resample(time='1D').mean()
daily_sst.plot()
plt.title('Daily Average SST at NDBC Buoy 44025')
plt.xlabel('')
plt.ylabel('Sea Surface Temperature (C)');

# Convert Dataset to Dataframe
df = ds.to_dataframe()
df = df.droplevel(['latitude','longitude']) # Drop extra indices that we don't need
df.head()

# Add a yearday column
df['yearday'] = df.index.dayofyear

# Calculate Annual Cycle
avg_sst = df.sea_surface_temperature.groupby(df.yearday).mean()

# Plot data by Yearday
plt.plot(df.yearday,df.sea_surface_temperature,'.',markersize=1,label='Raw Hourly Data');
avg_sst.plot(linewidth=3,label='10 Year Average')

plt.legend()
plt.xlabel('Day of Year')
plt.ylabel('Sea Surface Temperature (C)')
plt.title('Seasonal Cycle of Sea Surface Temperature at NDBC 44025 from 2010-2019');
plt.savefig('NDBC_44025_Seasonal_SST.png');

# Calculate daily average
df_daily = df.resample('1D').mean()
df_daily['yearday'] = df_daily.index.dayofyear
df_daily.head()

# Calculate SST Anomoly based on the 10-year average
df_daily['sst_climate'] = avg_sst[df_daily.yearday].values
df_daily['sst_anomoly'] = df_daily['sea_surface_temperature'] - df_daily['sst_climate']

# Anomoly Plot
df_daily['sea_surface_temperature'].plot(label='Raw Data')
df_daily['sst_climate'].plot(label='Climatic Prediction')
df_daily['sst_anomoly'].plot(label='Anomoly')
plt.legend(loc='upper left');
plt.xlabel('')
plt.ylabel('Sea Surface Temperature (C)')
plt.title('Sea Surface Temperature and Anomoly at NDBC 44025');

# Open the dataset and subset to the last 3 months
realtime = xr.open_dataset(url)
realtime = realtime.sel(time=slice('2019-01-01', '2020-04-01'))

# Convert to Dataframe
realtime = realtime.to_dataframe()
realtime = realtime.droplevel(['latitude','longitude']) # Drop extra indices

# Calculate daily average
realtime = realtime.resample('1D').mean()
realtime['yearday'] = realtime.index.dayofyear

# Add SST Anomoly based on the 10-year average
realtime['sst_climate'] = avg_sst[realtime.yearday].values
realtime['sst_anomoly'] = realtime['sea_surface_temperature'] - realtime['sst_climate']

# Finally, let's plot it!
fig = plt.figure(constrained_layout=True)

# Rearrange the subplots so we have 1 big and 1 small graph
gs = plt.GridSpec(nrows=3,ncols=1, figure=fig)
ax1 = fig.add_subplot(gs[0:2])
ax2 = fig.add_subplot(gs[2:3])

ax1.plot(realtime['sea_surface_temperature'],linewidth=2,label='Measured Data')
ax1.plot(realtime['sst_climate'],linewidth=2,label='Climatic Average')
ax2.plot(realtime['sst_anomoly'],linewidth=2,label='Anomoly')

ax1.legend(loc='upper left');
ax1.set_ylabel('Sea Surface Temperature (C)')
ax1.set_title('Sea Surface Temperature and Anomoly at NDBC 44025');
ax2.set_ylabel('Temperature (C)')
ax2.set_title('Anomoly')
plt.savefig('NDBC_44025_SST_Anomoly.png');