[Programming note: Office hours for all Bamboo Weekly paid subscribers will take place this coming Sunday. I’ll send e-mail about it in the coming 24 hours. Please come with your questions and comments!]
At first glance, the last week hasn’t seemed very encouraging when it comes to aviation safety. Between a crash in Japan and the side coming off of an Alaska Airlines plane, you might feel justified in thinking of airline travel as unsafe. In many ways, I would argue the opposite, that if no one was hurt in either of these accidents, we should be impressed by and grateful for the engineers who invest so much time improving airplane safety.
But what do the numbers say? Have there been fewer injuries over time or more? Do we see any specific trends?
This week, we looked at data collected by the National Transportation Safety Board (NTSB), part of the US Department of Transportation (https://www.ntsb.gov/Pages/home.aspx). They collect data about airline (and other transportation) accidents, and also investigate when an accident takes place. Their analysis helps recommend ways for the airline industry to improve and to avoid repeating mistakes.
Data and seven questions
This week's data comes from the NTSB's "CAROL" database (https://data.ntsb.gov/carol-main-public/basic-search), which lets you query all sorts of transportation incidents. You can then download the result of your query in either JSON or CSV format.
I asked you to download data from January 1st, 1998 through today from CAROL. I’ll admit that this required doing a bit of annoying work, since CAROL only lets you download up to 10,000 records at a time. Through a bit of experimentation, I found that I could download all of these records in JSON format by choosing “Aviation” as the search mode, specifying five pairs of dates, and then clicking on “search”. Then, after the search was over, I would click on “download JSON data.”
The zipfile that I got from each download contained a JSON-formatted file as well as a readme containing my query.
I chunked the data into these date ranges:
- 01/01/2018 to 12/31/2024
- 01/01/2012 to 12/31/2017
- 01/01/2007 to 12/31/2011
- 01/01/2002 to 12/31/2006
- 01/01/1998 to 12/31/2001
The files that I got were named:
- cases2024-01-10_08-14.json
- cases2024-01-10_08-15.json
- cases2024-01-10_08-17.json
- cases2024-01-10_08-18.json
- cases2024-01-10_08-20.json
From what I can tell, the filenames are based on the date when you downloaded the data, so you likely got files with different names.
This week, I gave you seven tasks and questions. Here are my answers and explanations; a link to the full Jupyter notebook is at the bottom of this message:
Combine the five JSON files into a single data frame. You will only need the following columns: `cm_mkey`, `cm_eventDate`, `cm_vehicles`, `cm_fatalInjuryCount`, `cm_seriousInjuryCount`, and `cm_minorInjuryCount`. Treat the `cm_eventDate` column as a date. Set the `cm_mkey` column to be the index. How can you be sure that you downloaded all of the files, covering all years?
First, and before doing anything else, I have to load Pandas:
import pandas as pdThen, with Pandas loaded, how can I turn the JSON files into a data frame? We’ve previously read numerous CSV files into Pandas using “read_csv”, and Excel files into Pandas using “read_excel”. The “read_json” function is similar, taking a filename as input and returning a data frame.
That’s fine for a single file and a single data frame, but how can we combine all five inputs into a single data frame? My preferred technique is to create a list of data frames using a list comprehension, invoking read_json on each file returned by “glob.glob”. I can then invoke “pd.concat” on that list of data frames, getting one data frame back:
import glob
all_dfs = [pd.read_json(one_filename)
for one_filename in glob.glob('cases*json')]
The good news is that this will indeed give me a list of data frames, one per file. However, I wanted to turn the “cm_eventDate” column into a datetime column. To do that, I pass the “convert_dates” keyword argument to read_json, passing a list of columns (one, in this case) that should be parsed in this way:
import glob
all_dfs = [pd.read_json(one_filename,
convert_dates=['cm_eventDate'])
for one_filename in glob.glob('cases*json')]
Now that we have a list of data frames, we can combine them with pd.concat:
df = (pd
.concat(all_dfs)
)The resulting data frame is quite wide, with 38 columns — most of which aren’t really necessary for our analysis. I thus decided to keep only those columns that I needed:
df = (pd
.concat(all_dfs)
[['cm_mkey', 'cm_eventDate', 'cm_vehicles',
'cm_fatalInjuryCount', 'cm_seriousInjuryCount', 'cm_minorInjuryCount']]
)I had asked that you set the “cm_mkey” column to be the index. Normally, I would have used the “index_col” keyword when reading the file from disk to choose an index column. But for reasons that I don’t quite understand, read_json doesn’t have such an option. I thus decided to handle it after creating the large data frame, invoking “set_index”:
df = (pd
.concat(all_dfs)
[['cm_mkey', 'cm_eventDate', 'cm_vehicles',
'cm_fatalInjuryCount', 'cm_seriousInjuryCount', 'cm_minorInjuryCount']]
.set_index('cm_mkey')
)I don’t know about you, but I made a number of mistakes when downloading the JSON files from the NTSB. How can we be sure that you haven’t forgotten any years when downloading the data?
I ran the following query:
(
df['cm_eventDate'].dt.year
.drop_duplicates()
.sort_values()
.diff()
)The above query grabs the year (thanks to the “dt.year” accessor on datetime columns) from each of the records, then invokes “drop_duplicates”, returning the unique years that are in the system. I sorted the years with “sort_values”, and calculated the difference between them with “diff”. The result was a series in which each row told me by how much it differed from the previous row.
My thinking was that if a year was missing, then I would have seen a gap of more than 1 between two years. But since I only had 1s in the result (except for NaN on the first row), I could be sure that I had data from all of the years.
Count the number of vehicles involved in each incident. How often does each count occur?
How often does an accident only involve one vehicle, and how often does it include more? If we had an integer column indicating how many vehicles were involved in an accident, then we could easily run “value_counts” and find out.
But instead, the “cm_vehicles” column contains Python objects — specifically, lists of dictionaries, with one dict per vehicle involved in each accident. We need to calculate the length of each list. Then we can take those numbers and run them through value_counts.
The solution is to use the “apply” method, which lets us run an arbitrary function on each element in a series. Here, we’ll run the “len” function, which will return the length of each list. We’ll get a series of integers back — one on which we can invoke value_counts. Here, I passed the normalize=True keyword argument, so that I would get percentages returned, rather than the raw integer values:
df['cm_vehicles'].apply(len).value_counts(normalize=True)The result:
cm_vehicles
1 0.985132
2 0.014805
3 0.000064
Name: proportion, dtype: float64I’ll add that this makes sense to me; 1-airplane accidents seem far more common than 2- or 3-airplane accidents, even though we know, thanks to the accident in Japan, that they do occur.
There is another way to do this, which is a bit sneaky — but which I’ve used on numerous occasions. Basically, the “str” accessor lets us run string methods on series that contain strings. So if “s” contains string data, I can call “s.str.len()” and get back a new series containing integers, each representing the string’s length.
The thing is, Pandas doesn’t check what underlying types are in a column of dtype “object”. If the underlying object supports the “len” function, then we can run it via “str.len()”, too. We can thus say:
df['cm_vehicles'].str.len().value_counts(normalize=True)I got the same results. Which raises the question of which is a better way to go. I used “%timeit” in my Jupyter notebook, and found that str.len ran in about 9 ms, whereas apply(len) ran in about 6.9 ms. So my sneaky trick of using str.len might be cool and interesting, but it’s also less efficient.
The `cm_vehicles` column contains a list of dictionaries. From the first element in each list, grab the `make`, `model`, and `operatorName` values and turn them into columns in the data frame.
As we’ve seen, the cm_vehicles column contains a list of dictionaries, with one dictionary for each vehicle in the accident. I asked you to add three columns to our data frame, one for each of the “make”, “model”, and “operatorName” keys in the first dictionary.
In other words, I want to apply something like “[0][‘make’]” to each of these lists, and assign the result to a new “make” column.
There are a few ways to do this, and in this case I decided to go with the sneaky trick that I mentioned earlier, namely using the “str” accessor. In particular, I decided to use “str.get”, which effectively invokes [] on a value. Here, I needed to use str.get twice — once for [0] to get the list’s first element, and again for the particular key.
I could have used str.get twice in a row. But instead, I decided to take advantage of “assign”, which lets us create new columns when we’re chaining methods. One of the nice things about assign is that it handles the keyword arguments in order. This means that later columns can refer to earlier columns, even in the same call to assign.
Let’s start by creating a new column, “first”, which will contain the first vehicle in the accident:
df = (
df
.assign(first=lambda df_: df_['cm_vehicles'].str.get(0))
)With that in place, we can now refer to the “first” column and create the three columns we really wanted to add. We’re still effectively using str.get twice in a row, but at least it’s a bit more readable:
df = (
df
.assign(first=lambda df_: df_['cm_vehicles'].str.get(0),
make=lambda df_: df_['first'].str.get('make'),
model=lambda df_: df_['first'].str.get('model'),
operatorName=lambda df_: df_['first'].str.get('operatorName'))
)Finally, since the “first” column is just there so that we can create the other three, I decided to remove it using “drop”, specifying that I want to drop a column, rather than a row:
df = (
df
.assign(first=lambda df_: df_['cm_vehicles'].str.get(0),
make=lambda df_: df_['first'].str.get('make'),
model=lambda df_: df_['first'].str.get('model'),
operatorName=lambda df_: df_['first'].str.get('operatorName'))
.drop('first', axis='columns')
)After running the above query, we now have three new columns with the make, model, and operator name of the first vehicle in the aviation accident.
What operator was involved in the greatest number of incidents each year?
Now that we have information about the operators — at least, the operator of the first vehicle in each incident — let’s find out which was involved in the greatest number of accidents per year.
Any time you hear “per year”, you know that we’ll want to perform a “groupby” operation. And any time you hear that we want to count on a per-value basis, you know that we want to run value_counts. Fortunately, value_counts is an aggregation method, so we can actually run it via groupby.
Moreover, we can run “groupby” on the years in our “cm_eventDate” column, thanks to the “dt.year” accessor:
(
df
.groupby(df['cm_eventDate'].dt.year)['operatorName']
.value_counts()
)Here’s what I got:
cm_eventDate operatorName
1998 1581
CONTINENTAL AIRLINES 4
TRANS WORLD AIRLINES 3
AMERICAN AIRLINES, INC. 3
UNKNOWN 2
...
2023 PABLO AIR CHARTERS LLC 1
2024 Japan Airlines 1
HANSEN JARED L 1
ALASKA AIRLINES INC 1
U.S Bangla Airlines 1
Name: count, Length: 19678, dtype: int64We have a series with a two-part multi-index. The outer part shows the years, and the inner part shows the number of incidents per operator in that year. That’s good, but we just want to know, per year, which operator had the maximum value.
One way to solve this is to turn our multi-indexed series into a pivot table. That is: We’ll move the years to be columns, and keep the inner part of the multi-index as the index. This is known as “unstacking,” and it’s a common Pandas idiom. We can even specify that we want to move the outer index level (aka, level=0) to be columns. The result:
(
df
.groupby(df['cm_eventDate'].dt.year)['operatorName']
.value_counts()
.unstack(level=0)
)We now have a data frame whose 17,489 rows represent all of the operators known to the database, and whose 27 columns represent the years in which there were incidents.
Now we can ask Pandas to find, for each year (column), the index that corresponds to the maximum value. We can do that with the “idxmax” method:
(
df
.groupby(df['cm_eventDate'].dt.year)['operatorName']
.value_counts()
.unstack(level=0)
.idxmax()
)The result:
cm_eventDate
1998
1999
2000
2001 AMERICAN AIRLINES INC
2002 Unknown
2003 Unknown
2004 On File
2005 On File
2006 On File
2007 On File
2008 Pilot
2009 Pilot
2010 Pilot
2011 Pilot
2012 Pilot
2013 Pilot
2014 Pilot
2015 Pilot
2016 Pilot
2017 Pilot
2018 Pilot
2019 Pilot
2020 Pilot
2021 Pilot
2022 Pilot
2023 UNITED AIRLINES INC
2024 ALASKA AIRLINES INC
dtype: objectFrom this, it would seem to me that yes, there have been incidents in which commercial airlines were responsible. For example, we can see the Alaska Airlines incident from last week already in the 2024 numbers. But for the most part, it seems like private, non-commercial pilots were far more likely to be involved in incidents than commercial airlines.
What 10 airplane makes + models were involved in the greatest number of incidents?
In the wake of the Alaska Airlines incident, there has been renewed talk about the safety of the Boeing 737 MAX. Have any particular makes and models of aircraft been involved in a high number of incidents?
I decided to find out by running “value_counts” on two of the columns, “make” and “model”:
df[['make', 'model']].value_counts().head(10)Here are the results I got:
make model
CESSNA 172 547
BOEING 737 501
Cessna 152 448
172 429
172N 320
172M 215
Piper PA-28-140 214
CESSNA 152 210
182 207
PIPER PA28 206
Name: count, dtype: int64I can see that the makes aren’t necessarily standardized for capitalization. Moreover, we can see that there are dashes and other characters in the model names that might not be standardized. I decided to use “str.lower” and “str.replace” on these columns, and then run value_counts again:
(
df[['make', 'model']]
.assign(make=lambda df_: df_['make'].str.lower(),
model=lambda df_: df['model'].str.replace('\W', '', regex=True))
.value_counts()
.head(10)
)Notice that I used a regular expression to search for any non-alphanumeric characters (\W), replacing them with empty strings. After running the above, here’s what I got:
make model
cessna 172 976
152 658
boeing 737 550
cessna 172N 512
piper PA28140 395
cessna 172S 386
182 377
piper PA18150 363
PA28 361
cessna 172M 358
Name: count, dtype: int64We can thus see that Cessna planes (models 172 and 152) were involved in far more incidents than others.
But wait a second: How much should we worry? (Not that I have a Cessna on back order, I’ll admit.) Maybe these are the most common makes and models to be in accidents, but it’s worth checking the percentages, thanks to normalize=True:
(
df[['make', 'model']]
.assign(make=lambda df_: df_['make'].str.lower(),
model=lambda df_: df['model'].str.replace('\W', '', regex=True))
.value_counts(normalize=True)
.head(10)
)The result:
make model
cessna 172 0.020735
152 0.013979
boeing 737 0.011685
cessna 172N 0.010877
piper PA28140 0.008392
cessna 172S 0.008201
182 0.008009
piper PA18150 0.007712
PA28 0.007669
cessna 172M 0.007606
Name: proportion, dtype: float64So yes, we see lots of Cessnas at the top of this list. But even the most common model to be involved in accidents, the Cessna 172, was involved in 2 percent of them.
Of course, we could just check makes rather than models:
df['make'].str.lower().value_counts(normalize=True).head(10)The results from here are a bit less encouraging for Cessna fans:
make
cessna 0.261222
piper 0.145857
beech 0.055209
boeing 0.042957
bell 0.029664
robinson 0.018877
mooney 0.014163
air tractor 0.009301
hughes 0.008515
bellanca 0.008366
Name: proportion, dtype: float64Wow. More than a quarter of all incidents in our data set involved Cessnas, and another 15 percent involved Pipers. That seems to fit the previous result we saw, which showed that the largest number of incidents in most years was listed as “pilot,” which I take to mean “private pilot.”
Has the number of injuries reported per year risen or declined over the years? Create a line plot from this data, showing all injury types.
I want to know whether the number of injuries reported per year has risen or fallen over time. I ran a “groupby” per year, asking for us to sum the count of fatal, serious, and minor injuries in each year:
(
df
.groupby(df['cm_eventDate'].dt.year)
[['cm_fatalInjuryCount', 'cm_seriousInjuryCount', 'cm_minorInjuryCount']].sum()
)The result is long, but a quick eyeballing of the data certainly looks like each type of injury has dropped. But it’ll definitely be easier to see it if we plot it:
(
df
.groupby(df['cm_eventDate'].dt.year)
[['cm_fatalInjuryCount', 'cm_seriousInjuryCount', 'cm_minorInjuryCount']].sum()
.plot.line()
)Here’s what I got:
It sure looks to me like there has been a massive decline in injuries over time, with a nice trend down. Of course, there have been upticks along the way, but for the most part, it looks like flying in 2024 is far safer than it was in 2000.
At what hour of the day do most flights with injuries (of any sort) occur?
Finally, I would guess that most flight problems and accidents take place at night. But is that true? Let’s find out how many injuries (total) occurred in each reported incident, and then find out how often each hour of the day has injuries.
I’ll start by creating a new column, total_injury_count, which sums all of the fatal, serious, and minor injuries. Notice that I’m using the “sum” method here, but because I want to sum across columns rather than down rows, I indicate that I want axis=“columns”:
(
df
.assign(total_injury_count = lambda df_: df_[['cm_fatalInjuryCount',
'cm_seriousInjuryCount',
'cm_minorInjuryCount']].sum(axis='columns'))
)With that in place, I can now perform a filter with “.loc”, keeping only those rows in which there was at least one injury:
(
df
.assign(total_injury_count = lambda df_: df_[['cm_fatalInjuryCount',
'cm_seriousInjuryCount',
'cm_minorInjuryCount']].sum(axis='columns'))
.loc[lambda df_: df_['total_injury_count'] > 0]
)Finally, I then grabbed the hour of the day with “dt.hour”, and counted how often each hour appeared:
(
df
.assign(total_injury_count = lambda df_: df_[['cm_fatalInjuryCount',
'cm_seriousInjuryCount',
'cm_minorInjuryCount']].sum(axis='columns'))
.loc[lambda df_: df_['total_injury_count'] > 0]
['cm_eventDate'].dt.hour.value_counts()
)The results:
cm_eventDate
14 1686
13 1672
16 1667
15 1647
12 1632
11 1625
17 1510
10 1433
18 1306
19 1201
9 1115
20 890
8 836
21 584
4 519
7 465
22 378
5 332
23 268
6 225
0 171
1 148
3 95
2 91
Name: count, dtype: int64Turns out, the most common hours for airline incidents is mid-day, from 10 a.m. - 5 p.m. I’m guessing (and this is just a guess) that this might have something to do with the fact that airport are busiest at these times of day, so the chances of an accident rise. But I’d be curious to hear thoughts and suggestions from you!
This week’s Jupyter notebook is at https://drive.google.com/file/d/1LkxmjdSiHN9Sxop7W3aPNx-thfOOiR-k/view?usp=sharing.
I’ll be back next Wednesday with more Python and Pandas puzzles based on current events.
Reuven