# data https://www.ncdc.noaa.gov/ibtracs/index.php?name=ib-v4-access
# ACE formula https://www.casact.org/pubs/forum/18spforumv2/05_Collins.pdf
# results validated against http://tropical.atmos.colostate.edu/Realtime/index.php?arch&loc=northatlantic
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


speedMin = [35, 55, 70, 85, 105, 125]
speedMax = [54, 69, 84, 104, 124, 999]
# Near-land trends since 1920 for correcting ACE
trends = [-0.00230, 0.000785, -0.00209, -0.001615, -0.001824, -0.001067]
colors = ['green', 'aqua', 'blue', 'purple', 'red', 'orange']
labels = ['35-54', '55-69', '70-84', '85-104', '105-124', '>=125']

YEARS = [1920, 2019]

def read_data(filename):
	ace = {}
	for speed in speedMin:
		ace[speed] = []
	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]
	aceYear = {}
	for speed in speedMin:
		aceYear[speed] = 0
	for line in csv_reader:
		firstToken = line[0]
		if firstToken.startswith('AL'):
			rows = int(line[columns1["rows"]])
			row = 0
			maxWind = 0
			aceStorm = 0
			isH = False
			isTS = False
			year = int(firstToken[4:9])
			if year > YEARS[1]:
				break
			if year < current_year:
				continue
			if year > current_year:
				for speed, trend in zip(speedMin, trends):
					# correct for detection/recording trend
					aceYear[speed] = aceYear[speed] * (1 - (trend * (YEARS[1] - current_year)))
					ace[speed].append(aceYear[speed])
				print(current_year, aceYear)
				for speed, trend in zip(speedMin, trends):
					aceYear[speed] = 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 are not in TS or HU status
			if stormType == 'TS' or stormType == 'HU':
				# ignore any times except 6 hour intervals,
				# ignore any intermediate reports with higher wind speeds
				if time == 0 or time == 600 or time == 1200 or time == 1800:
					aceStorm = aceStorm + (wind * wind)
					if wind > maxWind and stormType != 'EX':
						maxWind = wind
			row = row + 1
			if row == rows:
				# we are done with this storm
				for sMin, sMax in zip(speedMin, speedMax):
					if maxWind >= sMin and maxWind < sMax:
						aceYear[sMin] = aceYear[sMin] + (aceStorm / 10000)
	for speed, trend in zip(speedMin, trends):
		# correct for detection/recording trend
		aceYear[speed] = aceYear[speed] * (1 + (trend * (YEARS[1] - current_year)))
		ace[speed].append(aceYear[speed])
		aceYear[speed] = 0
	return ace
		
def plot():
	years = list(range(YEARS[0], YEARS[1] + 1))
	ace = read_data("hurdat.csv")

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

	aceTotal = []
	for year in years:
		aceYear = 0
		for speed, color, label in zip(speedMin, colors, labels):
			aceYear = aceYear + ace[speed][year - YEARS[0]]
		aceTotal.append(aceYear)

	plt.bar(years, aceTotal, color='red', label="ACE")

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

#	plt.tight_layout()
#	plt.legend(ncol = 4, loc='upper right')
	plt.title("ACE, corrected for detection trends")
	pngFile = "output.png"
	plt.savefig(pngFile)
	plt.show()

if __name__ == '__main__':	
	plot()