# https://www.ncdc.noaa.gov/ibtracs/index.php?name=ib-v4-access
import sys
import matplotlib.pyplot as plt
from matplotlib import gridspec
import numpy as np
from scipy import stats
import csv
from datetime import date

BASINTITLE = "ATLANTIC"
FILENAME = "hurdat.csv"
BASIN1 = "AL"
BASIN2 = "AL"
#BASINTITLE = "PACIFIC"
#FILENAME = "hurdat2-nepac.csv"
#BASIN1 = "EP"
#BASIN2 = "CP"

YEARS = [2010, 2019]

def read_data(filename):
	p35 = []
	p40 = []
	p45 = []
	p50 = []
	columns1 = {"desig": 0, "name": 1, "rows": 2}
	columns2 = {"ymd": 0, "time": 1, "landfall": 2, "stormType": 3, "lat": 4, "lon": 5, "wind": 6, "pressure": 7}
	f = open(filename, 'r')
	csv_reader = csv.reader(f, delimiter=',')
	current_year = YEARS[0]
	storms35 = 0
	storms40 = 0
	storms45 = 0
	storms50 = 0
	stormsTotal = 0
	for line in csv_reader:
		firstToken = line[0]
		if firstToken.startswith(BASIN1) or firstToken.startswith(BASIN2):
			rows = int(line[columns1["rows"]])
			row = 0
			maxWind = 0
			legitStorm = False
			year = int(firstToken[4:9])
			if year > YEARS[1]:
				break
			if year < current_year:
				continue
			if year > current_year:
				print(current_year, storms35, storms40, storms45, storms50, stormsTotal)
				p35.append(storms35 / stormsTotal)
				p40.append(storms40 / stormsTotal)
				p45.append(storms45 / stormsTotal)
				p50.append(storms50 / stormsTotal)
				storms35 = 0
				storms40 = 0
				storms45 = 0
				storms50 = 0
				stormsTotal = 0
				current_year = year
		elif year >= YEARS[0]:
			wind = int(line[columns2["wind"]])
			stormType = line[columns2["stormType"]].strip()
			# ignore storms that don't become TS or HU
			if stormType == 'TS' or stormType == 'HU':
				legitStorm = True
			if wind > maxWind and stormType != 'EX':
				maxWind = wind
			row = row + 1
			if row == rows:
				# we are done with this storm
				if legitStorm:
					stormsTotal = stormsTotal + 1
					if maxWind <= 50:
						storms50 = storms50 + 1
					if maxWind <= 45:
						storms45 = storms45 + 1
					if maxWind <= 40:
						storms40 = storms40 + 1
					if maxWind <= 35:
						storms35 = storms35 + 1
	print(current_year, storms35, storms40, storms45, storms50, stormsTotal)
	p35.append(storms35 / stormsTotal)
	p40.append(storms40 / stormsTotal)
	p45.append(storms45 / stormsTotal)
	p50.append(storms50 / stormsTotal)
	return p35, p40, p45, p50
		
def plot():
	years = list(range(YEARS[0], YEARS[1] + 1))
	p35, p40, p45, p50 = read_data(FILENAME)

	print(p35)

	fig = plt.figure(figsize=(8,4))
	plt.scatter(years, p35, color='green', label="peak35", alpha=0.6)
	plt.scatter(years, p40, color='blue', label="peak40", alpha=0.5)
	plt.scatter(years, p45, color='red', label="peak45", alpha=0.4)
	plt.scatter(years, p50, color='black', label="peak50", alpha=0.3)

	slope, intercept, r_value, p_value, std_err = stats.linregress(years, p35)
	points = [slope * y + intercept for y in years]
	plt.plot(years, points, color='green', alpha=0.5)

	slope, intercept, r_value, p_value, std_err = stats.linregress(years, p40)
	points = [slope * y + intercept for y in years]
	plt.plot(years, points, color='blue', alpha=0.5)
	
	slope, intercept, r_value, p_value, std_err = stats.linregress(years, p45)
	points = [slope * y + intercept for y in years]
	plt.plot(years, points, color='red', alpha=0.5)
	
	slope, intercept, r_value, p_value, std_err = stats.linregress(years, p50)
	points = [slope * y + intercept for y in years]
	plt.plot(years, points, color='black', alpha=0.5)

#	plt.tight_layout()
	plt.legend(ncol = 4, loc='upper right')
	plt.title(BASINTITLE + " percentage of storms by peak winds (knots) " + str(YEARS[0]) + "-" + str(YEARS[1]))
	pngFile = "output.png"
	plt.savefig(pngFile)
	plt.show()

if __name__ == '__main__':	
	plot()