#!/usr/bin/python

#start of week 2 work
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import sys, getopt

def main(argv):
	inputfile = ''
	location = ''
	try:
		opts, args = getopt.getopt(argv,"hi:o:",["ifile=","ofile="])
	except getopt.GetoptError:
		print('test.py -i <inputfile> -o <location>')
		sys.exit(2)
	for opt, arg in opts:
		if opt == '-h':
			print('test.py -i <inputfile> -o <location>')
			sys.exit()
		elif opt in ("-i", "--ifile"):
			inputfile = arg
		elif opt in ("-o", "--ofile"):
			location = arg
	# Use pandas to read in the csv file
	data_df = pd.read_csv("/Users/bruce/Desktop/shared/URI/OSC250/jupyter/" + inputfile)
	msl = data_df['MSL'] #mean sea level
	year = data_df['Year']
	plt.figure(figsize=(8,4))

	#scatterplot of year vs mean sea level
	plt.subplot(211)
	plt.scatter(year, msl, s=2, color='blue', alpha=0.5)
	plt.title('Time Series of Sea Level at ' + location)
	plt.xlabel('Year')
	plt.ylabel('Mean Sea Level')
	month = data_df['Month']

	#scatterplot of month vs mean sea level
	plt.subplot(212)
	plt.scatter(month, msl, s=1, color='green', alpha=0.5)
	time_x = year + month/12

	#scatterplot of detrended month vs mean sea level
	plt.ylabel('Mean Sea Level')
	plt.xlabel('Month')
	model = np.polyfit(time_x, msl, 1)
	predicted = np.polyval(model, time_x)
	detrended = msl - predicted

	plt.figure(figsize=(8,4))

	#scatterplot with trendline
	plt.subplot(131)
	plt.scatter(time_x, msl, s=1, color='#0000ff', alpha=0.5)
	plt.plot(time_x, predicted, linewidth=2, color="black")
	slope = model[0]
	intercept = model[1]
	plt.annotate('y$=%3.7sx+%3.7s$'%(slope, intercept), xy=(0.1, 0.9), xycoords='axes fraction', fontsize=8)
	plt.xlabel('Year', fontsize=14)
	r2=np.corrcoef(time_x, msl)[0,1]**2
	plt.annotate('$R^2=%3.4s$'%(r2), xy=(0.1,0.8),xycoords='axes fraction', fontsize=8)

	#detrended scatterplot
	plt.subplot(132)
	plt.scatter(time_x,detrended, s=1,color='red', alpha=0.5)
	plt.xlabel('Year', fontsize=14)
	plt.ylim(-2.0,2.0)

	#histogram
	plt.subplot(133)
	import matplotlib.mlab as mlab
	mu = np.mean(detrended)   #mean value
	sigma = np.std(detrended) #standard deviation
	x = mu + sigma * np.random.randn(5000) #artificial normal line for comparison
	n, bins, patches = plt.hist(x, 50, facecolor='green', alpha=0)
	plt.hist(detrended, 20, color='red', normed=1, orientation='horizontal')
	plt.xlabel('Count', fontsize=14)
	plt.ylim(-0.5,0.5)
	y = mlab.normpdf(bins, mu, sigma)
	l = plt.plot(y, bins, 'g--', linewidth=1)
	plt.tight_layout()
	plt.show()
	
if __name__ == "__main__":
   main(sys.argv[1:])

#https://tidesandcurrents.noaa.gov/waterlevels.html?id=8632200&units=standard&bdate=19200213&edate=20170209&timezone=GMT&datum=MSL&interval=m&action=data