This week, we looked at relief offered by FEMA, the Federal Emergency Management Agency (https://fema.gov), to a variety of disasters over the years. FEMA has provided food and water, as well as numerous other type of assistance, to people affected by Hurricane Helene, which passed through (among other places) the Southeastern United States.
I've heard about FEMA many times in the past, especially in the context of hurricanes, and I thought that it might be interesting to find out how often they're called upon to do their jobs, what sorts of disasters they deal with, and where the disasters typically occur.
Data and six questions
This week's data comes from FEMA, which like all US government agencies, publishes a great deal of data about its work. We can get a list of disasters on which it has worked, classified by type, date, and location (but not budget), from
https://www.fema.gov/openfema-data-page/disaster-declarations-summaries-v2
For this week's tasks, we'll use the Parquet data file that FEMA provides. Parquet is a binary format used by the Apache Arrow project (https://arrow.apache.org/). You can find the Parquet file (along with other formats) in the "Full data" section of the data page. You can download FEMA's Parquet file itself from:
https://www.fema.gov/api/open/v2/DisasterDeclarationsSummaries.parquet
We'll also use some data classified by FEMA regions. Information about those regions is at:
https://www.fema.gov/about/organization/regions
Here are my detailed solutions and explanations for this week. As always, a link to the Jupyter notebook I used to solve the problems is at the end of this message.
Here are my questions:
Create a data frame from the Parquet file. Make sure that all of the columns describing dates have a dtype of datetime.
We'll start by loading Pandas:
import pandas as pdWith that in place, we can then load the Parquet file:
filename = 'DisasterDeclarationsSummaries.parquet'
df = (
pd
.read_parquet(filename)
)It turns out that this works just fine; we get a data frame with all of the data in it. But upon closer examination, it seems that not all of the "date" columns were actually stored using datetime dtypes. That seems weird, but it makes sense if you consider that Parquet doesn't automatically set any types. Rather, it stores the data as it was told to store it. So if someone never made the effort to turn the "date" columns into datetime dtypes, Parquet will dutifully store them as strings.
To get around this, we'll need to overwrite the existing columns with datetime versions of themselves, invoking to_datetime on each of the columns.
The easiest way to do this, at least if we're using method chaining, is with the assign method. This lets us add a new column to the data frame on which we're working by passing a keyword argument. If the key matches an existing column, then we replace the column with the passed value – which can be a series or the result of a function call, such as a lambda.
In this case, I called assign with four keyword arguments. Each basically just replaces the existing column with a datetime version of the existing column:
filename = 'DisasterDeclarationsSummaries.parquet'
df = (
pd
.read_parquet(filename)
.assign(declarationDate =
lambda df_: pd.to_datetime(df_['declarationDate']),
incidentBeginDate =
lambda df_: pd.to_datetime(df_['incidentBeginDate']),
incidentEndDate =
lambda df_: pd.to_datetime(df_['incidentEndDate']),
disasterCloseoutDate =
lambda df_: pd.to_datetime(df_['disasterCloseoutDate']),
)
)With this in place, we have a data frame with 66,895 rows and 28 columns. All of the date-related columns now have the right dtype, too.
I should note that using a datetime rather than a string not only provides us with additional ability for calculations, but also reduces the size of the data, both on disk and in memory. The disk file we loaded is 6.3 MB in size, but after performing this transformation, it's only 5.7 MB in size. (I saved it with to_parquet after doing so.) And datetime values are only 64 bits, as opposed to strings, which are typically much larger. This is another example of how using an appropriate dtype isn't just nice, but really improves the performance of your system.
Create a bar graph showing the number of disasters that FEMA has handled in each year, starting in 2000. (Use the incidentBeginDate as the source of the year.) The bar should be subdivided into colored sections indicating the types of disasters. Only include the 5 most common disaster types. Does any year stand out? How different do things look if we count every entry in the data set as a separate disaster, vs. every unique disasterNumber?
I thought it would be interesting to see just how many disasters FEMA has to deal with each year, and to plot that, to see if the numbers have gone up or down over time. The best way to do this, I think, is to create a pivot table in which the index contains years and the columns contain disaster types, with the cells containing integers – the count of how many values match each type for each year.
Since FEMA deals with so many different types of disasters, I asked you to only look at those in the top 5. My favorite way to find the 5 most common elements is to run value_counts on the incidentType column, then get the top 5 values with head(5), then grab the index from there.
But I'm not interested in the top 5 incident types on their own. Rather, I only want those rows in the data frame whose incident type is in those top 5. I'll thus use isin inside of a .loc to keep only those rows whose incidentType is in that list:
(
df
.loc[lambda df_:
df_['incidentType'].isin(
df_['incidentType'].value_counts().head(5).index)]
)We'll need the year, and while we can get it with dt.year, that's not going to be good enough when we create a pivot table. So I decided to use assign to create a column called incidentBeginYear, containing just the year. And while I could then filter based on dt.year or the incidentBeginYear column, I figured that I had the column already, so why not use it? Here's how my query thus looks so far:
(
df
.loc[lambda df_:
df_['incidentType'].isin(
df_['incidentType'].value_counts().head(5).index)]
.assign(incidentBeginYear =
lambda df_: df_['incidentBeginDate'].dt.year)
.loc[lambda df_: df_['incidentBeginYear'] >= 2000]
)The data is now in place for us to create a pivot table. It's easiest to think of a pivot table as a 2-dimensional groupby, but working on two categorical columns (one for rows, the other for columns). We're interested in counting the number of intersections between these categories, so our aggfunc will be count, and we can choose any column we want for counting, so I'll just choose id. We then invoke pivot_table with the appropriate keyword arguments:
(
df
.loc[lambda df_:
df_['incidentType'].isin(
df_['incidentType'].value_counts().head(5).index)]
.assign(incidentBeginYear =
lambda df_: df_['incidentBeginDate'].dt.year)
.loc[lambda df_: df_['incidentBeginYear'] >= 2000]
.pivot_table(index='incidentBeginYear',
columns='incidentType',
values='id',
aggfunc='count')
)The resulting data frame shows us how many disasters of each type took place in each year. We can then create a bar plot with plot.bar. But that will give us a separate bar for each column of the data frame; I asked you to have a single bar for each year, divided into a colored segment for each incident type. To do that, we pass stacked=True as a keyword argument to plot.bar:
(
df
.loc[lambda df_:
df_['incidentType'].isin(
df_['incidentType'].value_counts().head(5).index)]
.assign(incidentBeginYear =
lambda df_: df_['incidentBeginDate'].dt.year)
.loc[lambda df_: df_['incidentBeginYear'] >= 2000]
.pivot_table(index='incidentBeginYear',
columns='incidentType',
values='id',
aggfunc='count')
.plot.bar(stacked=True)
)The result looks like this:
As you can see there was quite a rise in FEMA disasters in 2020. And whadaya know, the overwhelming majority were "Biological" in nature.
But note that this is counting each individual FEMA request as a separate disaster. That's not wrong, but this means that if a disaster is spread across multiple counties and states, we'll see numerous rows for each one. What if we just count each incident as a single FEMA event?
We can do that by invoking drop_duplicates on the data frame, indicating that we want to keep the disasterNumber column unique, keeping only the first row with each disaster ID number:
(
df
.drop_duplicates(subset='disasterNumber')
.loc[lambda df_:
df_['incidentType'].isin(
df_['incidentType'].value_counts().head(5).index)]
.assign(incidentBeginYear =
lambda df_: df_['incidentBeginDate'].dt.year)
.loc[lambda df_: df_['incidentBeginYear'] >= 2000]
.drop_duplicates(subset='disasterNumber')
.pivot_table(index='incidentBeginYear',
columns='incidentType',
values='id',
aggfunc='count')
.plot.bar(stacked=True)
)This results in a different graph:
According to this plot, 2020 certainly had a larger share of biological problems than most others, but it didn't dominate as much as in the other plot. Of course, perhaps the pandemic indicates that we should count each individual incident, because the notion that 2020 had fewer biological issues than 2011 seems, on the face of it, a bit ridiculous.
Create a histogram showing how frequently each FEMA disaster is declared, only counting each disasterNumber a single time. Are any months more (or less) likely to have disasters?
Next, I asked you to create a histogram of the months in which disasters took place. I started by again keeping each disaster number a single time:
(
df
.drop_duplicates(subset='disasterNumber')
)I then retrieved the month via the dt.month attribute:
(
df
.drop_duplicates(subset='disasterNumber')
['declarationDate'].dt.month
)Finally, I created a histogram with plot.hist, indicating that we should have 12 bins, to evenly look at the number of events per month:
(
df
.drop_duplicates(subset='disasterNumber')
['declarationDate'].dt.month
.plot.hist(bins=12)
)Here's the result:
We can see that the greatest number of incidents appear to take place in the summer and early fall, dropping off in October. I have to assume that this has a lot to do with hurricanes, severe storms, and tornadoes, all of which take place in the late summer and early fall. I'm not sure why March has a big climb in the number of incidents, but I'm open to suggestions.
We see that a single disaster number can be associated with a number of rows in the data set. By how much does that vary? What can we learn by comparing the mean and median?
We've now seen that each disaster number might be associated with a single row in the data set, or it might be associated with many. It all depends on how many separate locations asked for or were given assistance. I was curious to know by how much this varied.
We can find out by using groupby, counting the number of different values for id there are for each unique value of disasterNumber:
(
df
.groupby('disasterNumber')['id'].count()
)This returns a series whose index contains disaster numbers and whose values are the number of times that number appears in our data set. We can then run describe on the resulting series to get the descriptive statistics on the results:
(
df
.groupby('disasterNumber')['id'].count()
.describe()
)The result:
count 4981.000000
mean 13.430034
std 24.918811
min 1.000000
25% 1.000000
50% 3.000000
75% 14.000000
max 443.000000
Name: id, dtype: float64We can see that the number of rows per incident varies from 1 to 443. The mean is 13, but the median is 3 – and the 75% mark is only 14. This means that there is a small number of incident IDs with a very large number of rows. That's why it's so important to look at both the mean and median in such a data set; it can really tell you if the data is skewed up or down.
Which 10 states have had the greatest number of disasters (counting each record, not just each disaster number) starting in 2000?
Let's start by using dt.year , along with .loc and lambda, to only keep those rows in which the year is at least 2000:
(
df
.loc[lambda df_: df_['incidentBeginDate'].dt.year >= 2000]
)Next, we'll keep only the state column, then run value_counts on it, to find out how often each state appears. And then we'll use head to keep only the 10 most commonly mentioned states:
(
df
.loc[lambda df_: df_['incidentBeginDate'].dt.year >= 2000]
['state']
.value_counts()
.head(10)
)The result:
state
TX 3901
MO 2122
OK 2104
LA 2080
KY 1966
GA 1931
FL 1930
NC 1636
VA 1522
KS 1507
Name: count, dtype: int64We can see that Texas is far and away the leader in terms of FEMA-related disasters over the last quarter century, followed by Missouri and Oklahoma. Given all of the problems with fires and mudslides in California, I was sure that they could appear on this list, too.
FEMA divides the United States into 10 regions. Plot a bar graph showing the number of total disasters per region, starting in 2000, counting each disaster ID only once. The x axis should show a comma-separated list of states for each region, as found on the region-information page.
In the previous question, we looked at how often each state appeared in our FEMA data set. But FEMA divides the US into regions, and I was wondering if the regions were split more evenly. I thus asked you to create a bar plot showing how often each region appeared in the data set.
But there's a catch: I didn't want you to just show the region number. Rather, I wanted you to show the states in each region. That data is on the regional information page (https://www.fema.gov/about/organization/regions) that I mentioned in the "data" section, above.
To get that data, I used pd.read_html, a fantastic method that returns one data frame for each HTML table it finds at the URL we give:
(
pd.read_html('https://www.fema.gov/about/organization/regions')[0]
)read_html always returns a list of data frames. In this case, I know that we're using the data frame at index 0 on that list, which we grab as a data frame.
However, the data is a bit messy. The region contains the word "Region", and has some extra spaces in it, while we want integers. And the states are separated with | characters, when I asked for us to have comma-separated values.
Fortunately, we can use assign to correct this, along with some string-related methods from the str accessor. First, we can replace 'Region ' with the empty string, then calling astype to turn the result into an integer.
In the second case, we can use str.split to get back a Python list. Unlike the str.split that comes with Python, this one can use regular expressions – so we do that, looking for \s* (i.e., optional whitespace of any length), followed by | followed by more optional whitespace. We split on anything matching that to get our list.
We then immediately turn around and invoke str.join, putting ', ' (i.e., a comma and space) between the values).
This gives us everything we want, so we then drop the States / Territory column and set Region to be our index with set_index:
(
pd.read_html('https://www.fema.gov/about/organization/regions')[0]
.assign(Region = lambda df_:
df_['Region'].str.replace(
'Region ', '').astype(int),
states =
lambda df_: df_['States/ Territory'].str.split(
r'\s*\|\s*').str.join(', '))
.drop('States/ Territory', axis='columns')
.set_index('Region')
)Now that we have our regional data frame in place, we can join it with df. Note that join only works on the indexes of the two data frames, so we actually run join on df.set_index('region'). We then get them to match up, getting one large data frame in which the region ID number is the index, and the text describing states in that region is in the states column.
We then retrieve just the states column, giving us a series, and run value_counts on it. This tells us how many times each state appears. Finally, we create a horizontal bar plot with plot.barh:
(
pd.read_html('https://www.fema.gov/about/organization/regions')[0]
.assign(Region = lambda df_:
df_['Region'].str.replace(
'Region ', '').astype(int),
states =
lambda df_: df_['States/ Territory'].str.split(
r'\s*\|\s*').str.join(', '))
.drop('States/ Territory', axis='columns')
.set_index('Region')
.join(df.set_index('region'))
['states']
.value_counts()
.plot.barh()
)The result:
As you can see, the text was rather long, making the plot a bit weird looking. But no matter; we can now see which regions, and which states, get FEMA's attention most often. As it happens, the states with the greatest number of incidents are those that were affected by Hurricane Helene just over the last few days.
I hope that you enjoyed this! Here is a link to my Jupyter notebook: https://drive.google.com/file/d/1NUvt7Q32DVBTnJ-SFuK7G3XLyyFYbVZC/view?usp=sharing
I'll be back next Wednesday with more Pandas puzzles about current events and real-world data.
Reuven