# data https://www.ncdc.noaa.gov/ibtracs/index.php?name=ib-v4-access
# 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]
speedMax = [54, 69, 84]
# Near-land trends since 1920 for correcting count
trends = [-0.00230, 0.000785, -0.00209]
colors = ['green', 'aqua', 'blue']
labels = ['35-54', '55-69', '70-84']

#speedMin = [85, 105, 125]
#speedMax = [104, 124, 999]
# Near-land trends since 1920 for correcting count
#trends = [-0.001615, -0.001824, -0.001067]
#colors = ['purple', 'red', 'orange']
#labels = ['85-104', '105-124', '>=125']

YEARS = [1920, 2019]

def read_data(filename):
	count = {}
	for speed in speedMin:
		count[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]
	countYear = {}
	for speed in speedMin:
		countYear[speed] = 0
	for line in csv_reader:
		firstToken = line[0]
		if firstToken.startswith('AL'):
			rows = int(line[columns1["rows"]])
			row = 0
			maxWind = 0
			year = int(firstToken[4:9])
			if year > YEARS[1]:
				break
			if year < current_year:
				continue
			if year > current_year:
				print(current_year, countYear)
				for speed, trend in zip(speedMin, trends):
					# correct for detection/recording trend
					countYear[speed] = countYear[speed] * (1 - (trend * (YEARS[1] - current_year)))
					count[speed].append(countYear[speed])
				print(current_year, countYear)
				for speed in speedMin:
					countYear[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:
					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:
						countYear[sMin] = countYear[sMin] + 1
	for speed, trend in zip(speedMin, trends):
		# correct for detection/recording trend
		countYear[speed] = countYear[speed] * (1 + (trend * (YEARS[1] - current_year)))
		count[speed].append(countYear[speed])
	return count
		
def plot():
	years = list(range(YEARS[0], YEARS[1] + 1))
	count = read_data("hurdat.csv")

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

	for speed, color, label in zip(speedMin, colors, labels):
		plt.bar(years, count[speed], color=color, label=label, alpha=0.5)

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

#	plt.tight_layout()
	plt.legend(ncol = 3, loc='upper right')
	plt.title("hurricane count corrected for detection trend")
	pngFile = "output.png"
	plt.savefig(pngFile)
	plt.show()

if __name__ == '__main__':	
	plot()