Learn geopandas by plotting tornados on a map

In this tutorial we will take a look at the powerful geopandas library and use it to plot historical tornado data on a map of the United States.

We will use two different shapefiles from NOAA, the first dataset includes the origination point for each tornado. The second dataset includes a line path of each tornado. The data comes from the NOAA Storm Prediction Center an can be downloaded from https://www.spc.noaa.gov/gis/svrgis/

A quick note before we start

I assume you know some basic python and how to install jupyter to run the companion notebook. To get started, clone my git repository with the following commands.

git clone https://github.com/joncutrer/geopandas-tutorial
cd geopandas-tutorial
jupyter notebook

Here are the commands you will need to run if have not already install geopandas.

conda install geopandas
conda install -c conda-forge descartes

The data we will be working with comes from the US Census and is in a common shapefile format.

ls -al data

total 11044
drwxrwxr-x 2 sysadmin sysadmin    4096 May  3 19:45 ./
drwxrwxr-x 5 sysadmin sysadmin    4096 May  6 23:11 ../
-rw-rw-r-- 1 sysadmin sysadmin 3929143 Oct  1  2019 1950-2018-torn-aspath.zip
-rw-rw-r-- 1 sysadmin sysadmin 3191635 Oct  1  2019 1950-2018-torn-initpoint.zip
-rw-rw-r-- 1 sysadmin sysadmin 3808548 May  2  2019 cb_2018_48_tract_500k.zip
-rw-rw-r-- 1 sysadmin sysadmin       5 May  1 20:37 usa-states-census-2014.cpg
-rw-rw-r-- 1 sysadmin sysadmin   15201 May  1 20:37 usa-states-census-2014.dbf
-rw-rw-r-- 1 sysadmin sysadmin     143 May  1 20:37 usa-states-census-2014.prj
-rw-rw-r-- 1 sysadmin sysadmin     257 May  1 20:37 usa-states-census-2014.qpj
-rw-rw-r-- 1 sysadmin sysadmin  309672 May  1 20:37 usa-states-census-2014.shp
-rw-rw-r-- 1 sysadmin sysadmin   18517 May  1 20:37 usa-states-census-2014.shp.xml
-rw-rw-r-- 1 sysadmin sysadmin     564 May  1 20:37 usa-states-census-2014.shx

Import geopandas and load the U.S. map shapefile

import matplotlib.pyplot as plt
import geopandas
states = geopandas.read_file('data/usa-states-census-2014.shp')

states.head()

Reproject coordinates

By default this shp file is in WGS 84.

states.crs

Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

To make the map look a little more familiar lets reproject it’s coordinates to Mercator.

states = states.to_crs("EPSG:3395")

Now lets plot the U.S. map and see what we get

states.plot()

We can also plot the state polygons with no fill color by using boundary.

states.boundary.plot()

Load the NOAA US tornado dataset

Source: https://www.spc.noaa.gov/gis/svrgis/

One of the cool features of geopandas is that you can load a shapefile directly from a zip archive without first unzipping it. The file path to open a zipped shapefile is prefixed with zip://

However, there is a caveat. If the zip archive contains a subfolder you must specify it by appending !foldername to the path.

ls -al  data

total 11044
drwxrwxr-x 2 sysadmin sysadmin    4096 May  3 19:45 ./
drwxrwxr-x 5 sysadmin sysadmin    4096 May  6 23:11 ../
-rw-rw-r-- 1 sysadmin sysadmin 3929143 Oct  1  2019 1950-2018-torn-aspath.zip
-rw-rw-r-- 1 sysadmin sysadmin 3191635 Oct  1  2019 1950-2018-torn-initpoint.zip
-rw-rw-r-- 1 sysadmin sysadmin 3808548 May  2  2019 cb_2018_48_tract_500k.zip
-rw-rw-r-- 1 sysadmin sysadmin       5 May  1 20:37 usa-states-census-2014.cpg
-rw-rw-r-- 1 sysadmin sysadmin   15201 May  1 20:37 usa-states-census-2014.dbf
-rw-rw-r-- 1 sysadmin sysadmin     143 May  1 20:37 usa-states-census-2014.prj
-rw-rw-r-- 1 sysadmin sysadmin     257 May  1 20:37 usa-states-census-2014.qpj
-rw-rw-r-- 1 sysadmin sysadmin  309672 May  1 20:37 usa-states-census-2014.shp
-rw-rw-r-- 1 sysadmin sysadmin   18517 May  1 20:37 usa-states-census-2014.shp.xml
-rw-rw-r-- 1 sysadmin sysadmin     564 May  1 20:37 usa-states-census-2014.shx

Load the shapefile from zip archive

We will use the above mentioned technique to load the 1950-2018-torn-initpoint.shp shapefile from the zip archive.

tornados = geopandas.read_file(
    'zip://data/1950-2018-torn-initpoint.zip!1950-2018-torn-initpoint'
)

Note: If loading from inside a zip file confuses you don’t let is slow you down, just unzip the archive.

Lets have a look at the DataFrame

tornados.shape

(63645, 23)

tornados.info()

RangeIndex: 63645 entries, 0 to 63644
Data columns (total 23 columns):
 #   Column    Non-Null Count  Dtype   
---  ------    --------------  -----   
 0   om        63645 non-null  int64   
 1   yr        63645 non-null  int64   
 2   mo        63645 non-null  int64   
 3   dy        63645 non-null  int64   
 4   date      63645 non-null  object  
 5   time      63645 non-null  object  
 6   tz        63645 non-null  int64   
 7   st        63645 non-null  object  
 8   stf       63645 non-null  int64   
 9   stn       63645 non-null  int64   
 10  mag       63645 non-null  int64   
 11  inj       63645 non-null  int64   
 12  fat       63645 non-null  int64   
 13  loss      63645 non-null  float64 
 14  closs     63645 non-null  float64 
 15  slat      63645 non-null  float64 
 16  slon      63645 non-null  float64 
 17  elat      63645 non-null  float64 
 18  elon      63645 non-null  float64 
 19  len       63645 non-null  float64 
 20  wid       63645 non-null  int64   
 21  fc        63645 non-null  int64   
 22  geometry  63645 non-null  geometry
dtypes: float64(7), geometry(1), int64(12), object(3)
memory usage: 11.2+ MB

tornados

63645 rows × 23 columns

Reproject coordinates

We will also need to reproject it’s coordinates to Mercator so that our maps line up.

tornados = tornados.to_crs("EPSG:3395")

Take note how reprojecting our coordinates to Mercator moves the decimal point in our lat, long data. This will come into play later.

tornados[:1].geometry

0    POINT (-10043244.459 4662018.860)
Name: geometry, dtype: geometry

Our first plot of U.S. tornado data

tornados.plot( figsize=(12,9), color='red', marker='x', markersize=1)

Now lets plot the tornados data on top of the U.S. basemap¶

import matplotlib.pyplot as plt
fig = plt.figure(1, figsize=(12,9)) 
ax = fig.add_subplot()
states.boundary.plot(ax=ax, color='black', linewidth=.8)

tornados.plot(ax=ax, color='red', marker='.', markersize=1)

Wow, that’s great!

We even chose an ominous red v-shaped marker to represent the tornados. The two maps line up perfectly but it looks like we have some outliers data for Hawaii, Alaska, and Puerto Rico that are stretching out our map canvas.

We want to focus on the lower 48 states so now we need to narrow down our view. We can accomplish this in a couple of ways. One way is to remove any tornados from our dataframe that are not in the lower 48. Another approach is to simply limit the viewport area of the map with ax.set_xlim() and ax.set_ylim().

Our first thought might be to pass in typical lat,long data such as ax.set_xlim(-101.12345, 72.12345). After we reprojected our maps to Mercator our coordinates data now looks like this ax.set_xlim(-10112345, 7212345).

import matplotlib.pyplot as plt
fig = plt.figure(1, figsize=(12,9)) 
ax = fig.add_subplot()

# US Lower 48 Bounding Box
# -141.00000, 26.00000, -65.50000, 72.00000

ax.set_xlim(-14100244, -7200000)
ax.set_ylim(2600000, 6550000)

fig.suptitle('United States Tornado Map (1950-2018)', fontsize=16)

states.boundary.plot(ax=ax, color='black', linewidth=.8)

tornados.plot(ax=ax, color='red', marker='v', markersize=8)

Hey! That was a success, now let’s press on.

Get total tornados by State

I already know the answer to this question (because its my home) but… Which state do you think has had the most tornados?

Hint: In Texas we call them twisters

# Create a copy of the original DataFrame
twisters_by_state = tornados.copy()
# Add a new column and set the value to 1
twisters_by_state['tornados'] = 1

# use groupby() and count() to total up all the tornadoes by state
twisters_by_state = twisters_by_state[['st','tornados']].groupby('st').count()

# sort by most tornadoes first
twisters_by_state.sort_values('tornados', ascending=False).head(10)

Yep, Texas has the most tornados!

We can also easily plot this data by calling plot()

twisters_by_state.plot.bar(figsize=(12,6), title='Tornados by State (1950-2018)')

Notice that our sorting did not stick, this is because we did not reassign the DataFrame only output it to the screen.

We can chain in a call to .sort_values() to get our ordering to work in the plot.

I will also select the top 10 be adding [:10] which means “slice first 10 rows”.

Now we can see our state labels a little better.

twisters_by_state.sort_values('tornados', ascending=False)[:10].plot.bar(figsize=(12,6), title='Tornados by State (Top 10)')

If we just want a list of the top 10 states we can do something like this.

top10 = twisters_by_state.sort_values('tornados', ascending=False)[:10].index
top10

Index(['TX', 'KS', 'OK', 'FL', 'NE', 'IA', 'IL', 'MO', 'MS', 'CO'], dtype='object', name='st')

Plot tornados by state

Now that we know that Texas has had the most tornados lets subset our tornado dataset to only show Texas Twisters.

texas_map = states[states['NAME'] == 'Texas']
texas_map.plot()

texas_twisters = tornados[tornados['st'] == 'TX']
texas_twisters

8804 rows × 23 columns

texas_twisters.plot(markersize=1)

Now that we have both of our DataFrames texas_map and texas_twisters lets overlay combine and plot them.

import matplotlib.pyplot as plt
fig = plt.figure(1, figsize=(10,10)) 
ax = fig.add_subplot()

# Adjust bounding box to only Texas
# -119.00000, 29.30000, -103.80000, 44.00000

ax.set_xlim(-11910000, -10380000)
ax.set_ylim(2930000, 4400000)

fig.suptitle('Texas Tornados (1950-2018)', fontsize=16)

texas_map.boundary.plot(ax=ax, color='black', linewidth=.8)

texas_twisters.plot(ax=ax, color='red', marker='v', markersize=20)

Now lets move over to the path shapefile and see what we have to work with.

Load the NOAA US tornado paths dataset

tornado_paths = geopandas.read_file('zip://data/1950-2018-torn-aspath.zip!1950-2018-torn-aspath')

tornado_paths.head()

Reproject to Mercator

tornado_paths = tornado_paths.to_crs("EPSG:3395")

tornado_paths.plot(figsize=(14,8))

texas_twister_paths = tornado_paths[tornado_paths['st'] == 'TX']
texas_twister_paths.head()

5 rows × 23 columns

Plot tornado points and paths for Texas

Now that we have our tornado paths DataFrame narrowed down to Texas lets plot the paths.

For this next map I will plot the start point of each Tornado as pink and the path data as Red.

As you can see, path data does not exist for all recorded tornados.

import matplotlib.pyplot as plt
fig = plt.figure(1, figsize=(10,10)) 
ax = fig.add_subplot()

# Adjust bounding box to only Texas
# -119.00000, 29.30000, -103.80000, 44.00000

ax.set_xlim(-11910000, -10380000)
ax.set_ylim(2930000, 4400000)

fig.suptitle('Texas Tornado Paths (1950-2018)', fontsize=16)

texas_map.boundary.plot(ax=ax, color='black', linewidth=.8)
texas_twisters.plot(ax=ax, color='pink', marker='v', markersize=10)
texas_twister_paths.plot(ax=ax, color='red')

Color coded tornados by year

In this variation of the Texas tornado map we are color coding the tornado markers and paths by year with a colormap file. To do this we add the parameter column='yr' to the call to plot().

The colormap used cmap='coolwarm' allows us to see older tornado records as cool blue colors and more recent tornados as warm red colors. Note that we also removed the color parameter as cmap and color cannot be used simultaneously.

import matplotlib.pyplot as plt
fig = plt.figure(1, figsize=(10,10)) 
ax = fig.add_subplot()

# Adjust bounding box to only Texas
# -119.00000, 29.30000, -103.80000, 44.00000

ax.set_xlim(-11910000, -10380000)
ax.set_ylim(2930000, 4400000)

fig.suptitle('Texas Tornado Paths (1950-2018)', fontsize=16)

texas_map.boundary.plot(ax=ax, color='black', linewidth=.8)
texas_twisters.plot(ax=ax, column='yr', cmap="coolwarm", marker='v', markersize=10)
texas_twister_paths.plot(ax=ax, column='yr', cmap="coolwarm")

Congratulations, you are well on your way to becoming a GIS expert with geopandas. I hope you have enjoyed this Python Tutorial, let me know by leaving a comment below.

If you enjoyed this article you might also like Learn geopandas by plotting a United States shapefile