# 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
from elevations import Elevations
from math import sin, cos, sqrt, atan2, radians

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

TOP10WETTEST = [{'rank': 1, 'total': 60.58, 'name': 'HARVEY', 'year': 2017, 'location': 'Nederland, TX'},
				{'rank': 2, 'total': 48.00, 'name': 'AMELIA', 'year': 1978, 'location': 'Medina, TX'},
				{'rank': 3, 'total': 45.20, 'name': 'EASY', 'year': 1950, 'location': 'Yankeetown, FL'},
				{'rank': 4, 'total': 45.00, 'name': 'CLAUDETTE', 'year': 1979, 'location': 'Alvin, TX'},
				{'rank': 5, 'total': 43.31, 'name': 'IMELDA', 'year': 2019, 'location': 'Jefferson County, TX'},
				{'rank': 6, 'total': 40.68, 'name': 'ALLISON', 'year': 2001, 'location': 'NW Jefferson County, TX'},
				{'rank': 7, 'total': 38.46, 'name': 'GEORGES', 'year': 1998, 'location': 'Munson, FL'},
				{'rank': 8, 'total': 36.71, 'name': 'DANNY', 'year': 1997, 'location': 'Dauphin Island Sea Lab, AL'},
				{'rank': 9, 'total': 35.93, 'name': 'FLORENCE', 'year': 2018, 'location': 'Elizabethtown, NC'},
				{'rank': 10, 'total': 29.76, 'name': 'UNNAMED-01', 'year': 1960, 'location': 'Port Lavaca #2, TX'}]

META = [{'year': 1950, 'name': 'EASY', '24h': ' ', 'total': 45.2, 'location': 'FL', 'notes': 'state record'},
		{'year': 1960, 'name': 'UNNAMED-01', '24h': 22, 'total': ' ', 'location': ' ', 'notes': 'KY state record (11.25)'},
		{'year': 1962, 'name': 'TD3', '24h': 22, 'total': ' ', 'location': 'LA', 'notes': 'state 24h record'},
		{'year': 1963, 'name': 'FLORA', '24h': 28.39, 'total': 100.39, 'location': 'Cuba', 'notes': '2nd highest in western hemisphere'},
		{'year': 1972, 'name': 'AGNES', '24h': '', 'total': ' ', 'location': '', 'notes': 'PA TC record (19.0)'},
		{'year': 1979, 'name': 'CLAUDETTE', '24h': 42, 'total': ' ', 'location': '', 'notes': 'US 24h record'},
		{'year': 1979, 'name': 'FREDERIC', '24h': '', 'total': ' ', 'location': '', 'notes': 'OH TC record (8.67)'},
		{'year': 1980, 'name': 'HERMINE', '24h': '', 'total': 31.15, 'location': 'Mexico', 'notes': ''},
		{'year': 1982, 'name': 'CHRIS', '24h': '', 'total': ' ', 'location': '', 'notes': 'TN TC record (13.6)'},
		{'year': 1994, 'name': 'ALBERTO', '24h': ' ', 'total': 27.85, 'location': 'GA', 'notes': 'state record'},
		{'year': 1997, 'name': 'DANNY', '24h': 32.52, 'total': 36.71, 'location': 'AL', 'notes': 'state record'},
		{'year': 1998, 'name': 'MITCH', '24h': ' ', 'total': 62.87, 'location': 'Honduras, Nicaragua', 'notes': '11,000 fatalities'},
		{'year': 1998, 'name': 'GEORGES', '24h': ' ', 'total': 38.46, 'location': 'FL/MS', 'notes': 'MS state record (32.21)'},
		{'year': 2000, 'name': 'KEITH', '24h': ' ', 'total': 32.67, 'location': 'Belize', 'notes': 'record'},
		{'year': 2001, 'name': 'ALLISON', '24h': ' ', 'total': '', 'location': '', 'notes': 'LA 2nd highest'},
		{'year': 2005, 'name': 'DENNIS', '24h': ' ', 'total': 42.99, 'location': 'Cuba', 'notes': '2nd highest'},
		{'year': 2005, 'name': 'WILMA', '24h': 64.33, 'total': ' ', 'location': '', 'notes': 'record 24h Western Hemisphere'},
		{'year': 2007, 'name': 'NOEL', '24h': ' ', 'total': 29.43, 'location': 'Bahamas', 'notes': 'record'},
		{'year': 2017, 'name': 'HARVEY', '24h': ' ', 'total': 60.58, 'location': 'TX', 'notes': 'TX and US record'},
		{'year': 2018, 'name': 'FLORENCE', '24h': ' ', 'total': 35.93, 'location': 'NC/SC', 'notes': 'state records'},
		{'year': 2019, 'name': 'BARRY', '24h': ' ', 'total': ' ', 'location': '', 'notes': 'AR state record (16.59)'},
		{'year': 2019, 'name': 'DORIAN', '24h': ' ', 'total': 22.84, 'location': 'Bahamas', 'notes': '2nd highest'}]

# NOTE: all labels below need to be unique text since they are used as keys
STORMS = {'min': [35, 55, 70, 85, 105, 125],
		  'max': [54, 69, 84, 104, 124, 999],
		  # Near-land trends since 1920 for correcting counts and ACE
		  'trends': [-0.00230, 0.000785, -0.00209, -0.001615, -0.001824, -0.001067],
		  'colors': ['green', 'aqua', 'blue', 'purple', 'red', 'black'],
		  'labels': ['35-54', '55-69', '70-84', '85-104', '105-124', '>=125']}

RI24 = {'min': [30, 40, 50, 60],
		'colors': ['blue', 'purple', 'red', 'black'],
		'labels': ['>=30', '>=40', '>=50', '>=60']}

STARTINGWINDS = {'min': [0, 35, 50, 65, 80, 105],
				 'max': [34, 49, 64, 79, 104, 999],
				 'colors': ['orange', 'lightgray', 'darkgray', 'gray', 'dimgray', 'black'],
				 'labels': ['<35', '35-49', '50-64', '65-79', '80-104', '>=105']}

SLOWING = {'title': "Percent of landfalling storms with slow motion ",
		   'max': [3, 6, 9, 12],
		   'colors': ['black', 'purple', 'blue', 'green'],
		   'labels': ['0-3', '0-6', '0-9', '0-12']}

elevations = Elevations();
elevations.loadFiles();

# HURDAT2 has two types of rows, first for each storm, second for each storm report
COLUMNS1 = {"desig": 0, "name": 1, "rows": 2}
COLUMNS2 = {"ymd": 0, "time": 1, "landfall": 2, "stormType": 3, "lat": 4, "lon": 5, "wind": 6, "pressure": 7}

# slow moving landfall table
class SpeedTable():
	def __init__(self, lower, upper):
		self.lowerThresh  = lower
		self.upperThresh = upper
		self.table = "<TABLE id=\"HURTABLE\"><THEAD><TR><TH>Year</TH><TH>Name</TH><TH>speed<br>mph</TH><TH>rank<br>US</TH><TH>Rain<br>24 hr</TH><TH>Rain<br>total</TH><TH>Notes</TH></TR></THEAD>"

	def addYear(self, year):
		stormCount = 0
		for stormNum in year.storms:
			storm = year.storms[stormNum]
			if storm['slowestSpeedNearLand'] <= self.upperThresh:
				stormCount = stormCount + 1
		if stormCount == 0:
			return
		firstRow = True
		for stormNum in year.storms:
			storm = year.storms[stormNum]
			if storm['slowestSpeedNearLand'] <= self.upperThresh:
				if firstRow:
					self.table = self.table + "<TBODY>\n\t<TR><TD rowspan=" + str(stormCount) + ">" + str(year.year) + "</TD>"
					firstRow = False
				else:
					self.table = self.table + "\t<TR>"
				if storm['slowestSpeedNearLand'] <= self.lowerThresh:
					self.table = self.table + "<TD class=\"highlight\">"
				else:
					self.table = self.table + "<TD>"
				#Year Name speed USrank 24hrRain totalRain Notes
				self.table = self.table + storm['name'] + "</TD><TD>" + str(storm['slowestSpeedNearLand']) + "</TD>"
				foundTop10 = False
				for top10 in TOP10WETTEST:
					if top10['year'] == year.year and top10['name'] == storm['name']:
						foundTop10 = True
						break
				if foundTop10:
					self.table = self.table + "<TD>" + str(top10['rank']) + "</TD>"
				else:
					self.table = self.table + "<TD></TD>"
				foundMeta = False
				for row in META:
					if row['year'] == year.year and row['name'] == storm['name']:
						foundMeta = True
						break
				if foundMeta:
					self.table = self.table + "<TD>" + str(row['24h']) + "</TD>"
					self.table = self.table + "<TD>" + str(row['total']) + "</TD>"
					self.table = self.table + "<TD>" + str(row['location']) + ' ' + str(row['notes']) + "</TD>"
				else:
					self.table = self.table + "<TD></TD><TD></TD><TD></TD>"
				self.table = self.table + "</TR>\n"
		self.table = self.table + "</TBODY>\n"

	def fini(self):
		self.table = self.table + "</TABLE>"
		return self.table

class Year():
	def __init__(self, yearNum):
		self.year = yearNum
		self.storms = {}
		self.stormNum = 0

	def addStorm(self, storm):
		self.storms[self.stormNum] = storm.summary
		self.stormNum = self.stormNum + 1

# FROM STORMS:
#{'desig': 'AL141990', 'madeTS': True, 'madeHU': True, 'firstWind': 45, 'lastWind': 40, 'maxW': 65, 'ace': 5.7875, 'n24HU': 6, 'n24RI': 0, 'max6Inc': 5, 'max24Inc': 15, 'max36Inc': 15, 'max6Dec': -10, 'max24Dec': -15, 'windAtLandfall': 0, 'wind6AfterLandfall': 0, 'wind24AfterLandfall': 0, 'slowestSpeedNearLand': 999}
#{'desig': 'AL151990', 'madeTS': True, 'madeHU': False, 'firstWind': 25, 'lastWind': 15, 'maxW': 55, 'ace': 1.4475, 'n24HU': 0, 'n24RI': 0, 'max6Inc': 5, 'max24Inc': 20, 'max36Inc': 25, 'max6Dec': -10, 'max24Dec': -20, 'windAtLandfall': 25, 'wind6AfterLandfall': 0, 'wind24AfterLandfall': 0, 'slowestSpeedNearLand': 6.990032819711232}
#{'desig': 'AL161990', 'madeTS': True, 'madeHU': True, 'firstWind': 25, 'lastWind': 20, 'maxW': 75, 'ace': 6.1725, 'n24HU': 6, 'n24RI': 0, 'max6Inc': 10, 'max24Inc': 30, 'max36Inc': 40, 'max6Dec': -10, 'max24Dec': -35, 'windAtLandfall': 0, 'wind6AfterLandfall': 0, 'wind24AfterLandfall': 0, 'slowestSpeedNearLand': 999}
# TO YEAR:
#1990 {'numTS': 14, 'numHU': 8, 'total24HU': 76, 'total24RI': 0, 'ACE': 91.1675, '35-54': 3, '55-69': 4, '70-84': 4, '85-104': 2, '105-124': 1, '>=125': 0, '>=30': 4, '>=40': 0, '>=50': 0, '>=60': 0, '<35': 15, '35-49': 1, '50-64': 0, '65-79': 0, '80-104': 0, '>=105': 0, '0-3': 0, '0-6': 0, '0-9': 0, '0-12': 0}
	def summarize(self):
		self.summary = {'numTS': 0, 'numHU': 0, 'numLF': 0, 'total24HU': 0, 'total24RI': 0, 'ACE': 0}
		for label in STORMS['labels']:
			self.summary[label] = 0
		for label in RI24['labels']:
			self.summary[label] = 0
		for label in STARTINGWINDS['labels']:
			self.summary[label] = 0
		for label in SLOWING['labels']:
			self.summary[label] = 0
		for stormNum in self.storms:
			storm = self.storms[stormNum]
			if storm['madeTS']:
				self.summary['numTS'] = self.summary['numTS'] + 1
			if storm['madeHU']:
				self.summary['numHU'] = self.summary['numHU'] + 1
			if storm['windAtLandfall'] > 0:
				self.summary['numLF'] = self.summary['numLF'] + 1
			self.summary['ACE'] = self.summary['ACE'] + storm['ace']
			self.summary['total24HU'] = self.summary['total24HU'] + storm['n24HU']
			self.summary['total24RI'] = self.summary['total24RI'] + storm['n24RI']
			for windMin, windMax, label in zip(STORMS['min'], STORMS['max'], STORMS['labels']):
				if storm['maxW'] >= windMin and storm['maxW'] <= windMax:
					self.summary[label] = self.summary[label] + 1
			for riMin, label in zip(RI24['min'], RI24['labels']):
				if storm['max24Inc'] >= riMin:
					self.summary[label] = self.summary[label] + 1
			for swMin, swMax, label in zip(STARTINGWINDS['min'], STARTINGWINDS['max'], STARTINGWINDS['labels']):
				if storm['firstWind'] >= swMin and storm['firstWind'] <= swMax:
					self.summary[label] = self.summary[label] + 1
			for slowMax, label in zip(SLOWING['max'], SLOWING['labels']):
				if storm['slowestSpeedNearLand'] <= slowMax:
					self.summary[label] = self.summary[label] + 1
		#print(self.year, self.summary)

class Storm():
	def __init__(self, line):
		self.desig = line[COLUMNS1["desig"]]
		self.name = line[COLUMNS1["name"]].strip()
		if self.name == "UNNAMED":
			self.name = self.name + '-' + self.desig[2:4]
		self.nrows = int(line[COLUMNS1["rows"]])
		self.rows = {}
		self.rowNum = 0
		
	# https://stackoverflow.com/questions/19412462/getting-distance-between-two-points-based-on-latitude-longitude
	def getDistance(self, lat1, lon1, lat2, lon2):
		# approximate radius of earth in miles
		R = 3960
		lat1 = radians(lat1)
		lon1 = radians(lon1)
		lat2 = radians(lat2)
		lon2 = radians(lon2)
		dlon = lon2 - lon1
		dlat = lat2 - lat1
		a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
		c = 2 * atan2(sqrt(a), sqrt(1 - a))
		distance = R * c
		return distance

	def convertLatitude(self, latString):
		if latString[-1] == 'N':
			return float(latString[0:-1])
		else:
			return float(latString[0:-1]) * -1.0
		
	def convertLongitude(self, lonString):
		if lonString[-1] == 'W':
			return float(lonString[0:-1]) * -1.0
		else:
			return float(lonString[0:-1])

	def addStormRow(self, line):
		time = int(line[COLUMNS2["time"]])
		if time not in [0, 600, 1200, 1800]:
			return
		lat = self.convertLatitude(line[COLUMNS2["lat"]])
		lon = self.convertLongitude(line[COLUMNS2["lon"]])
		ele = elevations.getElevation(lon, lat)
		overland = False
		if ele != None and ele > -500:
			overland = True
		wind = int(line[COLUMNS2["wind"]])
		ymd = int(line[COLUMNS2["ymd"]])
		stormType = line[COLUMNS2["stormType"]].strip()
		self.rows[self.rowNum] = {'lat': lat, 'lon': lon, 'overland': overland, 'wind': wind, 'type': stormType}
		self.rowNum = self.rowNum + 1

#{'lat': 19.1, 'lon': -85.7, 'overland': False, 'wind': 45, 'type': 'TS'}
#{'lat': 19.7, 'lon': -85.5, 'overland': False, 'wind': 50, 'type': 'TS'}
#{'lat': 20.2, 'lon': -85.4, 'overland': False, 'wind': 60, 'type': 'TS'}
#{'lat': 20.9, 'lon': -85.1, 'overland': False, 'wind': 65, 'type': 'HU'}
#{'lat': 21.7, 'lon': -85.1, 'overland': False, 'wind': 75, 'type': 'HU'}
#{'lat': 22.7, 'lon': -85.2, 'overland': False, 'wind': 85, 'type': 'HU'}
#{'lat': 23.7, 'lon': -85.8, 'overland': False, 'wind': 85, 'type': 'HU'}
#{'lat': 24.6, 'lon': -86.2, 'overland': False, 'wind': 90, 'type': 'HU'}
#{'lat': 25.6, 'lon': -86.4, 'overland': False, 'wind': 100, 'type': 'HU'}
#{'lat': 26.6, 'lon': -86.5, 'overland': False, 'wind': 110, 'type': 'HU'}
#{'lat': 27.7, 'lon': -86.6, 'overland': False, 'wind': 120, 'type': 'HU'}
#{'lat': 29.0, 'lon': -86.3, 'overland': False, 'wind': 125, 'type': 'HU'}
#{'lat': 30.2, 'lon': -85.4, 'overland': True, 'wind': 135, 'type': 'HU'}
#{'lat': 31.5, 'lon': -84.5, 'overland': True, 'wind': 80, 'type': 'HU'}
#{'lat': 32.8, 'lon': -83.2, 'overland': True, 'wind': 50, 'type': 'TS'}
#{'lat': 34.1, 'lon': -81.7, 'overland': True, 'wind': 45, 'type': 'TS'}
#{'lat': 35.6, 'lon': -80.0, 'overland': True, 'wind': 45, 'type': 'TS'}
#{'lat': 36.5, 'lon': -77.7, 'overland': True, 'wind': 50, 'type': 'EX'}
	def summarizeStorm(self):
		self.summary = {'desig': self.desig, 'name': self.name, 'madeTS': False, 'madeHU': False,
						'firstWind': 0, 'lastWind': 0, 'maxW': 0, 'ace': 0,
						'n24HU': 0, 'n24RI': 0,
						'max6Inc': 0, 'max24Inc': 0, 'max36Inc': 0,
						'max6Dec': 0, 'max24Dec': 0,
						'windAtLandfall': 0, 'wind6AfterLandfall': 0, 'wind24AfterLandfall': 0,
						'slowestSpeedNearLand': 999}

		# point to the last 6 rows 
		agos = [36, 30, 24, 18, 12, 6]
		rowsAgo = {}
		for ago in agos:
			rowsAgo[ago] = None

		for rownum in self.rows:
			row = self.rows[rownum]
			#print(row)
			if self.summary['firstWind'] == 0:
				self.summary['firstWind'] = row['wind']
			if row['wind'] > self.summary['maxW']:
				self.summary['maxW'] = row['wind']
			if row['type'] == 'HU':
				self.summary['madeHU'] = True
			if row['type'] == 'TS':
				self.summary['madeTS'] = True

			# Official ACE: only use 6 hour intervals with TS or HU
			# validated with https://weather.fandom.com/wiki/1992_Atlantic_hurricane_season
			if row['type'] == 'HU' or row['type'] == 'TS':
				self.summary['ace'] = self.summary['ace'] + (row['wind'] * row['wind'])

			if row['type'] != 'EX':
				# Rapid Intensification
				if rowsAgo[6] != None and row['wind'] - rowsAgo[6]['wind'] > self.summary['max6Inc']:
					self.summary['max6Inc'] = row['wind'] - rowsAgo[6]['wind']
				if rowsAgo[24] != None and row['wind'] - rowsAgo[24]['wind'] > self.summary['max24Inc']:
					self.summary['max24Inc'] = row['wind'] - rowsAgo[24]['wind']
				if rowsAgo[36] != None and row['wind'] - rowsAgo[36]['wind'] > self.summary['max36Inc']:
					self.summary['max36Inc'] = row['wind'] - rowsAgo[36]['wind']

				# Rapid Weakening, excluding over land
				if rowsAgo[6] != None and row['overland'] == False and rowsAgo[6]['overland'] == False and row['wind'] - rowsAgo[6]['wind'] < self.summary['max6Dec']:
					self.summary['max6Dec'] = row['wind'] - rowsAgo[6]['wind']
				if rowsAgo[24] != None and row['overland'] == False and rowsAgo[6]['overland'] == False and rowsAgo[12]['overland'] == False and rowsAgo[18]['overland'] == False and rowsAgo[24]['overland'] == False and row['wind'] - rowsAgo[24]['wind'] < self.summary['max24Dec']:
					self.summary['max24Dec'] = row['wind'] - rowsAgo[24]['wind']

				if self.summary['windAtLandfall'] == 0 and row['overland'] == True:
					self.summary['windAtLandfall'] = row['wind']

				# allow one TS reading before HU periods as a compromise for RI ratio
				if rowsAgo[24] != None and (rowsAgo[24]['type'] == 'HU' or rowsAgo[24]['type'] == 'TS') and rowsAgo[18]['type'] == 'HU' and rowsAgo[12]['type'] == 'HU' and rowsAgo[6]['type'] == 'HU' and row['type'] == 'HU':
					self.summary['n24HU'] = self.summary['n24HU'] + 1
					if row['wind'] - rowsAgo[24]['wind'] > 30:
						self.summary['n24RI'] = self.summary['n24RI'] + 1

				# find the slowest 6 hour movement near landfall, but leave out TD only
				if self.summary['madeTS']:
					previousLat = 0
					previousLon = 0
					minDistance = 999
					nearLand = False
					for ago in agos:
						if rowsAgo[ago] != None:
							lat = rowsAgo[ago]['lat']
							lon = rowsAgo[ago]['lon']
							if rowsAgo[ago]['overland']:
								nearLand = True
							distance = self.getDistance(previousLat, previousLon, lat, lon)
							if distance < minDistance:
								minDistance = distance
							previousLat = lat
							previousLon = lon
					if nearLand and minDistance < 999:
						if self.summary['slowestSpeedNearLand'] > minDistance / 6:
							self.summary['slowestSpeedNearLand'] = round(minDistance / 6, 2)

			# shift the pointers to the last six rows and point to the current row
			for ago in agos:
				if ago > 6:
					rowsAgo[ago] = rowsAgo[ago - 6]
			rowsAgo[6] = row
		self.summary['ace'] = self.summary['ace'] / 10000
		self.summary['lastWind'] = row['wind']
		#print(self.summary)

def read_data(filename, startYear, endYear):
	years = []
	f = open(filename, 'r')
	csv_reader = csv.reader(f, delimiter=',')
	current_year = startYear
	year = Year(current_year)
	for line in csv_reader:
		firstToken = line[0]
		if firstToken.startswith(BASIN1) or firstToken.startswith(BASIN2):
			rows = int(line[COLUMNS1["rows"]])
			row = 0
			storm = Storm(line)
			yearNum = int(firstToken[4:9])
			# done with years we want
			if yearNum > endYear:
				break
			# skip over years we don't want
			if yearNum < current_year:
				continue
			# finished the current year
			if yearNum > current_year:
				year.summarize()
				years.append(year)
				current_year = yearNum
				year = Year(current_year)
		elif yearNum == current_year:
			storm.addStormRow(line)
			row = row + 1
			if row == rows:
				storm.summarizeStorm()
				year.addStorm(storm)
	year.summarize()
	years.append(year)
	return years
		
def plotSpeed():
	yearStart = 1950
	yearEnd = 2019
	years = read_data(FILENAME, yearStart, yearEnd)

	table = SpeedTable(3, 6)
	for year in years:
		table.addYear(year)
	html = table.fini()
	f = open("slowstorms.html", "w")
	f.write(html)
	f.close()

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

	for slowMax, color, label in zip(SLOWING['max'], SLOWING['colors'], SLOWING['labels']):
		validYears = []
		validPs = []
		for year in years:
			p = year.summary[label] / year.summary['numLF'] * 100
			validYears.append(year.year)
			validPs.append(p)

		axes = plt.gca()
		plt.scatter(validYears, validPs, color=color, label=label, alpha=0.6)

		if len(validYears) > 0:
			slope, intercept, r_value, p_value, std_err = stats.linregress(validYears, validPs)
			print(slowMax, slope)
			points = [slope * y + intercept for y in validYears]
			plt.plot(validYears, points, color=color, alpha=0.5)

	plt.legend(ncol = 4, loc='upper right')
	plt.title(SLOWING['title'] + str(yearStart) + "-" + str(yearEnd))
	pngFile = "slow" + str(yearStart) + ".png"
	plt.savefig(pngFile)
	plt.show()

if __name__ == '__main__':	
	plotSpeed()