Virus Surge
Interactive map with range selector (slider) showing the spread of virus hotspots across the USA over time.
Today we will make this brilliant visualization in the NYT Article U.S. Virus Cases Climb Toward a Third Peak. Both the static visualization and the interactive range slider.
import pandas as pd
import geopandas as gpd
import altair as alt
alt.renderers.set_embed_options(actions=False)
Shapefiles for US States from the topojson/us-atlas repository. Shapefiles from US Census are fine too.
states_shp_uri = 'https://cdn.jsdelivr.net/npm/us-atlas@3/states-10m.json'
states = gpd.read_file(states_shp_uri)
states['lon'] = states['geometry'].centroid.x
states['lat'] = states['geometry'].centroid.y
We will remove the states that are not supported in the albersUsa projection.
# Albers USA projection does not support the following.
states = states[~((states['id'] == '69') | (states['id'] == '78') | (states['id'] == '60') | (states['id'] == '72') | (states['id'] == '66'))]
states.head()
Initial Visualization
alt.Chart(states).mark_geoshape().encode().project('albersUsa')
Getting the COVID data as usual from NYT GitHub Repository
us_covid_data_uri = 'https://github.com/nytimes/covid-19-data/blob/master/us-counties.csv?raw=true'
raw_data = pd.read_csv(us_covid_data_uri)
raw_data['date'] = pd.to_datetime(raw_data['date'])
covid = raw_data.copy()
covid.head()
We will also load the shapefile for counties because we need them for calculating the longitude and latitudes of the centroids of polygons of the counties
counties_shp_uri = 'https://cdn.jsdelivr.net/npm/us-atlas@3/counties-10m.json'
counties = gpd.read_file(counties_shp_uri)
counties.head()
Removing any empty geometries beacause that will give errors when calculating longitudes and latitudes from the centroid
counties.geometry.is_empty.any()
empty = counties.geometry.is_empty
counties = counties[~empty]
counties['lon'] = counties['geometry'].centroid.x
counties['lat'] = counties['geometry'].centroid.y
counties.head()
counties.rename(columns={'id': 'fips'}, inplace=True)
counties.fips = counties.fips.astype(int)
The usual NYC aggregation because of the way data is reported by NYT.
# Extract the boroughs in shapefile geodataframe for NYC
boroughs_nyc = counties[counties['fips'].isin(['36005', '36047', '36061', '36081', '36085'])]
boroughs_nyc['State'] = "NYC"
#combined_nyc = boroughs_nyc.dissolve(by='STATEFP')
agg_nyc_data = boroughs_nyc.dissolve(by='State').reset_index()
agg_nyc_data['fips'] = 1
agg_nyc_data['lon'] = agg_nyc_data['geometry'].centroid.x
agg_nyc_data['lat'] = agg_nyc_data['geometry'].centroid.y
agg_nyc_data = agg_nyc_data.drop(columns=['State'])
counties = gpd.GeoDataFrame(pd.concat([counties, agg_nyc_data], ignore_index=True))
counties.tail()
covid = raw_data.copy()
covid = covid[~((covid['state'] == "Puerto Rico")|(covid['state'] == "Virgin Islands")|(covid['state'] == "Guam")|(covid['state'] == "Northern Mariana Islands"))]
# Extract the boroughs in shapefile geodataframe for NYC
boroughs_nyc = counties[counties['fips'].isin(['36005', '36047', '36061', '36081', '36085'])]
boroughs_nyc['State'] = "NYC"
#combined_nyc = boroughs_nyc.dissolve(by='STATEFP')
agg_nyc_data = boroughs_nyc.dissolve(by='State').reset_index()
agg_nyc_data['fips'] = 1
agg_nyc_data['lon'] = agg_nyc_data['geometry'].centroid.x
agg_nyc_data['lat'] = agg_nyc_data['geometry'].centroid.y
agg_nyc_data = agg_nyc_data.drop(columns=['State'])
counties = gpd.GeoDataFrame(pd.concat([counties, agg_nyc_data], ignore_index=True))
# Make fips in covid data for "New York City" as 1 to reflect what we have done above
covid.loc[covid['county'] == 'New York City','fips'] = 1
#covid = covid.merge(counties[['fips', 'lon', 'lat']], on='fips', how='inner')
#covid = covid.drop(columns=['STATEFP', 'geometry'])
# Getting per day statistics on CASES
covid = covid.assign(daily=covid.groupby('fips')['cases'].diff())
#covid = covid[covid['date'] > '2020-05-01'] # Working with data from May
# Getting past 2 weeks data for cases per given date - sum of cases in past 2 weeks
covid = covid.assign(past_2_weeks = covid.groupby('fips')['daily'].apply(lambda case: case.rolling(window=14).sum().shift(1)))
# Keeping only the latest date data
covid = covid[covid['date'] == '2020-10-27']
# Merging population estimate data - we could have done it earlier but what the heck!
## Getting census data
census = pd.read_csv('co-est2019-alldata.csv')
census = census[['STATE', 'COUNTY', 'STNAME', 'CTYNAME','POPESTIMATE2019']] # Keeping only 2019 population estimates
## Constructing FIPS from STATE and COUNTY codes in census itself
census['STATE'] = census.STATE.astype(str)
census['COUNTY'] = census.COUNTY.astype(str)
census['STATE'] = census.STATE.str.pad(2, fillchar='0')
census['COUNTY'] = census.COUNTY.str.pad(3, fillchar='0')
census = census.assign(fips = census.STATE + census.COUNTY)
census.fips = census.fips.astype(int)
census = census.drop(columns=['STATE','COUNTY','STNAME','CTYNAME'])
## Merging the census data with covid data on FIPS
covid = covid.merge(census, on='fips', how='inner')
## Getting the per capita data
covid = covid.assign(past_2_weeks_per_capita = (covid.past_2_weeks/covid.POPESTIMATE2019)*1000)
covid = covid[covid['past_2_weeks_per_capita'].notna()]
# Per Capita Deaths
covid = covid.assign(deaths_per_capita = covid.deaths/covid.POPESTIMATE2019*1000)
# Per Capita Deaths
covid = covid.assign(cases_per_capita = covid.cases/covid.POPESTIMATE2019*1000)
# Finally merging with geometry
plot_data = counties.merge(covid, on='fips', how='inner')
#plot_data = covid
#plot_data = plot_data.drop(columns=['STATEFP', 'daily', 'past_2_weeks', 'past_2_weeks_per_capita', 'deaths_per_capita'])
print(plot_data.info())
plot_data.head()
Making a custom height column with SVG path strings so that we scale the data using that column
plot_data = plot_data.assign(height = plot_data['past_2_weeks_per_capita'].apply(lambda x: f"M -1 0 L 0 -{x*2} L 1 0"))
Plotting the data -
#collapse
base = alt.Chart(states).mark_geoshape(fill='#ededed', stroke='white').encode(text='name:N', longitude='lon:Q', latitude='lat:Q').project('albersUsa')
text = alt.Chart(states).mark_text().encode(text='name:N', longitude='lon:Q', latitude='lat:Q', tooltip=['name:N']).project('albersUsa')
spikes = alt.Chart(plot_data).mark_point(
fillOpacity=0.5,
fill='red',
strokeOpacity=1,
strokeWidth=1,
stroke='red',
).encode(
shape=alt.Shape('height', scale=None),
longitude='lon:Q',
latitude='lat:Q',
).project('albersUsa').properties(width=1200, height=900)
(base+text+spikes).configure_view(stroke=None)