| Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli |
This file and the rest of the files in this directory were copied from https://github.com/applied-math/demos and adapted slightly to run on ScienceData.
Preparations¶
In [3]:
# This will take a long while, please have patience...
import os
os.system("conda install cartopy")
Out[3]:
In [4]:
%matplotlib inline
import numpy
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
In [19]:
import warnings
import cartopy
warnings.filterwarnings("ignore")
Computational Mathematics¶
Study of how we use computers to solve equations. This includes:
- Predicting the weather and climate change
- Design of cars, airplanes and rockets
- Data science and machine learning
All of these applications require the use of
- Linear algebra
- Statistics
- Calculus
Studying Hurricanes¶
Today we are going to do 3 things:
- Compute a model of the wind and pressure around a hurricane
- Visualize the tracks and strengths of past storms in the Atlantic
- Compute some statistics of hurricanes in the Atlantic
Today we're going to visualize hurricanes by find an equation that modelsusing historical tracks of hurricanes in the Atlantic ocean since 1950. Take a look at our dataset.
Hurricanes and Typhoons Around the World¶
Models of Wind and Pressure in a Hurricane¶
One way we can guess at what the wind and pressure around a hurricane. One of the ways we can do that is to use the following equations.
Wind $$ W(r) = \sqrt{100 \frac{A \cdot B}{\rho_{\text{air}}} \frac{e^{-A / r^B}}{r^B} + \frac{r^2 \cdot f^2}{4}} - \frac{r \cdot f}{2} $$
Pressure $$ P(r) = P_c + (P_n - P_c) e^{-A / r^B} $$
In [5]:
# Parameters
N = 1000
radius = 100e3 # Radius of storm
Pn = 1005 # Normal pressure
Pc = 950 # Central pressure
A = 23.0 # Shape of the storm
B = 1.5
rho_air = 1.15 # Density of air
OMEGA = 7.2722052166430395e-5 # Rotation of the earth (constant)
theta_0 = 0.52359877559829882 # Latitude (in radians) that we are at
f = 2.0 * OMEGA * numpy.sin(theta_0) # Coriolis parameter
# Evaluate profiles
x = numpy.concatenate((numpy.linspace(-radius,-0.01,N), numpy.linspace(0.01,radius,N)),axis=0)
r = numpy.abs(x) * 1e-3
p = Pc + (Pn - Pc) * numpy.exp(-A/(r)**B)
C = 1e1**2 * A * B * (Pn-Pc) / rho_air
v = numpy.sqrt(C * numpy.exp(-A / r**B) / r**B + r**2 * f**2 / 4.0) - r * f / 2.0
x /=1e3
fig = plt.figure()
fig.set_figwidth(fig.get_figwidth() * 2)
axes = fig.add_subplot(1, 2, 1)
axes.plot(x, v, color='blue')
axes.set_title("Wind Velocity Profile")
axes.set_xlabel('km')
axes.set_ylabel('m/s')
axes.set_xlim([numpy.min(x), numpy.max(x)])
axes.set_ylim([0.0, numpy.max(v) + 5])
axes = fig.add_subplot(1, 2, 2)
axes.plot(x, p, color='blue')
axes.set_title("Pressure Profile")
axes.set_xlabel('km')
axes.set_ylabel('mb')
axes.set_xlim([numpy.min(x), numpy.max(x)])
axes.set_ylim([Pc - 5, Pn + 5])
plt.show()
Visualization of Historical Hurricanes¶
In [20]:
#simple map
plt.figure(figsize=(15,9))
central_lat = 30
central_lon = -70
ax = plt.axes(projection=ccrs.LambertConformal(central_lon, central_lat))
ax.set_extent([-100, -40, 5, 57])
ax.add_feature(cfeature.LAND, edgecolor='black')
gl=ax.gridlines(draw_labels = True, x_inline=False, y_inline=False)
gl.top_labels = False
ax.set_title('Atlantic Hurricane Tracks (1950-2012)')
plt.show()
Play around with the code. Here some questions to help you get started.
- What happens when you delete
mapping.drawcoastlines()from the map? - What happens when you change the figure size?
- Can you change the colors?
Loading Data¶
First we need to load the data from the text file.
In [22]:
hurricanes = {}
os.system("curl -LO https://sciencedata.dk/public/6e3ed434c0fa43df906ce2b6d1ba9fc6/analyzing_hurricanes/hurdat2-1851-2020-052921.txt")
with open("./hurdat2-1851-2020-052921.txt", 'r') as hurdat_data:
for line in hurdat_data:
split_data = line.split(",")
# Header line
if split_data[0][:2] == "AL":
name = split_data[1].lstrip()
N = int(split_data[2])
if name == "UNNAMED":
# Skip this storm
for i in range(N):
line = hurdat_data.readline()
else:
i = 0
while name in hurricanes.keys():
i += 1
name += str(i)
if name in hurricanes.keys():
name = name[:-1]
hurricanes[name] = {'year': None,
'month': numpy.empty(N, dtype=int),
'day': numpy.empty(N, dtype=int),
'hour': numpy.empty(N, dtype=int),
'latitude': numpy.empty(N),
'longitude': numpy.empty(N),
'wind_speed': numpy.empty(N),
'pressure': numpy.empty(N),
'category': numpy.empty(N, dtype=int)}
# Read data
for i in range(N):
line = hurdat_data.readline()
if len(line.strip()) == 0:
break
split_data = line.split(",")
hurricanes[name]['year'] = int(split_data[0][:4])
hurricanes[name]["month"][i] = int(split_data[0][4:6])
hurricanes[name]["day"][i] = int(split_data[0][6:8])
hurricanes[name]["hour"][i] = int(split_data[1])
hurricanes[name]["latitude"][i] = float(split_data[4][:-1])
hurricanes[name]["longitude"][i] = -float(split_data[5][:-1])
hurricanes[name]["wind_speed"][i] = float(split_data[6])
hurricanes[name]["pressure"][i] = float(split_data[7])
# Categorize hurricanes (wind speed in knots)
if (hurricanes[name]["wind_speed"][i] >= 64 and
hurricanes[name]["wind_speed"][i] < 83):
hurricanes[name]["category"][i] = 1
elif (hurricanes[name]["wind_speed"][i] >= 83 and
hurricanes[name]["wind_speed"][i] < 96):
hurricanes[name]["category"][i] = 2
elif (hurricanes[name]["wind_speed"][i] >= 96 and
hurricanes[name]["wind_speed"][i] < 113):
hurricanes[name]["category"][i] = 3
elif (hurricanes[name]["wind_speed"][i] >= 113 and
hurricanes[name]["wind_speed"][i] < 135):
hurricanes[name]["category"][i] = 4
elif hurricanes[name]["wind_speed"][i] >= 135:
hurricanes[name]["category"][i] = 5
else:
hurricanes[name]["category"][i] = 0
# Count hurricanes
base_year = 3000
for name in hurricanes.keys():
base_year = min(base_year, hurricanes[name]['year'])
num_hurricanes = numpy.zeros(2021 - base_year, dtype=int)
for name in hurricanes.keys():
num_hurricanes[hurricanes[name]['year'] - base_year] += 1
Now let's plot the data onto the map.¶
In [23]:
# Plot the data
fig = plt.figure(figsize=(15,9))
central_lat = 30
central_lon = -70
ax = plt.axes(projection=ccrs.LambertConformal(central_lon, central_lat))
ax.set_extent([-100, -40, 5, 57])
ax.add_feature(cfeature.LAND, edgecolor='black')
gl=ax.gridlines(draw_labels = True, x_inline=False, y_inline=False)
gl.top_labels = False
ax.set_title('Atlantic Hurricane Tracks (1950-2021)')
for (name, data) in iter(hurricanes.items()):
longitude = data['longitude']
latitude = data['latitude']
ax.plot(longitude, latitude, linewidth=1.5, color='r', transform = ccrs.Geodetic())
plt.show()
In [25]:
# Plot the data
category_color = {5:'red', 4:'yellow', 3:'orange', 2:'green', 1:'blue', 0:'gray'}
plt.figure(figsize=(15,9))
central_lat = 30
central_lon = -70
ax = plt.axes(projection=ccrs.LambertConformal(central_lon, central_lat))
ax.set_extent([-100, -40, 5, 57])
ax.add_feature(cfeature.LAND, edgecolor='black')
gl=ax.gridlines(draw_labels = True, x_inline=False, y_inline=False)
gl.top_labels = False
ax.set_title('Atlantic Hurricane Tracks (1950-2021)')
for (name, data) in iter(hurricanes.items()):
longitude = data['longitude']
latitude = data['latitude']
for i in range(len(longitude) - 1):
color = category_color[data['category'][i]]
ax.plot(longitude[i:i+2], latitude[i:i+2], linewidth=1.5, color=color, transform = ccrs.Geodetic())
for (category, color) in iter(category_color.items()):
ax.plot([0], [0], color=color, label="Category %s" % category)
ax.legend(loc=2)
plt.show()
The data is still pretty hard to discern as it looks like colorful spaghetti. We can ask more specific questions such as
- Which storms happened in 2010?
- What storms were category 3 or above?
- Where did Hurricane Sandy go?
In [26]:
# Plot the data
category_color = {5:'red', 4:'yellow', 3:'orange', 2:'green', 1:'blue', 0:'gray'}
plt.figure(figsize=(15,9))
central_lat = 30
central_lon = -70
ax = plt.axes(projection=ccrs.LambertConformal(central_lon, central_lat))
ax.set_extent([-100, -40, 5, 57])
ax.add_feature(cfeature.LAND, edgecolor='black')
gl=ax.gridlines(draw_labels = True, x_inline=False, y_inline=False)
gl.top_labels = False
ax.set_title('Atlantic Hurricane Tracks (1950-2012)')
for (name, data) in iter(hurricanes.items()):
if data['year'] == 2011 and numpy.any(numpy.array(data['category']) >= 3):
longitude = data['longitude']
latitude = data['latitude']
for i in range(len(longitude) - 1):
color = category_color[data['category'][i]]
ax.plot(longitude[i:i+2], latitude[i:i+2], linewidth=1.5, color=color, transform = ccrs.Geodetic())
for (category, color) in iter(category_color.items()):
ax.plot([0], [0], color=color, label="Category %s" % category)
ax.legend(loc=2)
plt.show()
Statistics of Hurricanes¶
Now let us look at some statistics related to the number of hurricanes in the Atlantic each year from 1950 onward. We will also look at the mean and a "windowed" mean of the data.
In [27]:
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)
axes.plot(range(1950, 2021, 1), num_hurricanes, 'o-')
axes.plot((1950, 2013), (num_hurricanes.mean(), num_hurricanes.mean()), 'r')
# Compute rolling average
num_hurricanes_mean = numpy.zeros(15)
for i in range(0, 2021-1950, 5):
num_hurricanes_mean[i//5] = num_hurricanes[i:i+5].mean()
axes.plot(range(1950, 2021, 5), num_hurricanes_mean, 'ko--')
axes.set_xlim((1950, 2021))
axes.set_title("Hurricanes Each Year")
axes.set_xlabel("Year")
axes.set_ylabel("Number")
plt.show()
Recap¶
- What did we learn to do today?
- What are some issues that came up?
- Below are some more datasets. Take a look and see what kinds of questions you can come up that you want to explore?
The Collection of Really Great, Interesting, Situated Datasets