# 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 = "
Year | Name | speed mph | rank US | Rain 24 hr | Rain total | Notes |
"
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 + "\n\t" + str(year.year) + " | "
firstRow = False
else:
self.table = self.table + "\t
"
if storm['slowestSpeedNearLand'] <= self.lowerThresh:
self.table = self.table + ""
else:
self.table = self.table + " | "
#Year Name speed USrank 24hrRain totalRain Notes
self.table = self.table + storm['name'] + " | " + str(storm['slowestSpeedNearLand']) + " | "
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 + "" + str(top10['rank']) + " | "
else:
self.table = self.table + " | "
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 + "" + str(row['24h']) + " | "
self.table = self.table + "" + str(row['total']) + " | "
self.table = self.table + "" + str(row['location']) + ' ' + str(row['notes']) + " | "
else:
self.table = self.table + " | | | "
self.table = self.table + "
\n"
self.table = self.table + "\n"
def fini(self):
self.table = self.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()