# 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

YEARS = [2000, 2019]

def read_data(filename):
	pct30 = []
	pct40 = []
	pct50 = []
	pct60 = []
	hurricanes = []
	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]
	num30Year = 0
	num40Year = 0
	num50Year = 0
	num60Year = 0
	hurrYear = 0
	for line in csv_reader:
		firstToken = line[0]
		if firstToken.startswith('AL'):
			rows = int(line[columns1["rows"]])
			row = 0
			max24Storm = 0
			priorWinds = [0, 0, 0, 0]
			priorYMD = [0, 0, 0, 0]
			priorTimes = [9999, 9999, 9999, 9999]
			isHurricane = False
			year = int(firstToken[4:9])
			if year > YEARS[1]:
				break
			if year < current_year:
				continue
			if year > current_year:
				print(current_year, num30Year, num40Year, num50Year, num60Year, hurrYear)
				if hurrYear > 0:
					pct30.append(num30Year / hurrYear * 100)
					pct40.append(num40Year / hurrYear * 100)
					pct50.append(num50Year / hurrYear * 100)
					pct60.append(num60Year / hurrYear * 100)
				else:
					pct30.append(0)
					pct40.append(0)
					pct50.append(0)
					pct60.append(0)
				hurricanes.append(hurrYear)
				num30Year = 0
				num40Year = 0
				num50Year = 0
				num60Year = 0
				hurrYear = 0
				current_year = year
		elif year >= YEARS[0]:
			wind = int(line[columns2["wind"]])
			ymd = int(line[columns2["ymd"]])
			time = int(line[columns2["time"]])
			stormType = line[columns2["stormType"]].strip()
			# ignore storms that don't become HU
			if stormType == 'HU':
				isHurricane = True
			if stormType == 'TS' or stormType == 'HU':
				if time == 0 or time == 600 or time == 1200 or time == 1800:
					if priorWinds[0] > 0 and time == priorTimes[0] and ymd == priorYMD[0] + 1:
						wind24 = wind - priorWinds[0]
						if wind24 > max24Storm:
							max24Storm = wind24
					priorWinds = priorWinds[1:] + [wind]
					priorYMD = priorYMD[1:] + [ymd]
					priorTimes = priorTimes[1:] + [time]
			row = row + 1
			if row == rows:
				# we are done with this storm
				if isHurricane:
					hurrYear = hurrYear + 1
					if max24Storm >= 30:
						num30Year = num30Year + 1
					if max24Storm >= 40:
						num40Year = num40Year + 1
					if max24Storm >= 50:
						num50Year = num50Year + 1
					if max24Storm >= 60:
						num60Year = num60Year + 1
	if hurrYear > 0:
		pct30.append(num30Year / hurrYear * 100)
		pct40.append(num40Year / hurrYear * 100)
		pct50.append(num50Year / hurrYear * 100)
		pct60.append(num60Year / hurrYear * 100)
	else:
		pct30.append(0)
		pct40.append(0)
		pct50.append(0)
		pct60.append(0)
	hurricanes.append(hurrYear)
	return hurricanes, pct30, pct40, pct50, pct60
		
def plot():
	years = list(range(YEARS[0], YEARS[1] + 1))
	hurricanes, pct30, pct40, pct50, pct60 = read_data("hurdat.csv")

	fig = plt.figure(figsize=(8,4))
	plt.scatter(years, pct30, color='yellow', label=">= 30 kt")
	plt.scatter(years, pct40, color='green', label=">= 40 kt")
	plt.scatter(years, pct50, color='blue', label=">= 50 kt")
	plt.scatter(years, pct60, color='red', label=">= 60 kt")
	plt.bar(years, hurricanes, width=1.0, alpha=0.5, color='gray')

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

	slope, intercept, r_value, p_value, std_err = stats.linregress(years, pct40)
	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, pct50)
	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, pct60)
	points = [slope * y + intercept for y in years]
	plt.plot(years, points, color='red', alpha=0.5)

#	plt.tight_layout()
	plt.legend(ncol = 4, loc='upper right')
	plt.title("percent of storms with >= 30/40/50/60 kt 24h wind increase")
	pngFile = "output.png"
	plt.savefig(pngFile)
	plt.show()

if __name__ == '__main__':	
	plot()