Week 4A: More Interactive Data Viz, Working with Raster Data¶
Sep 27, 2021
Housekeeping¶
- Homework #2 due on Wednesday (9/29)
- Choose a dataset to visualize and explore
- OpenDataPhilly or one your choosing
- Email me if you want to analyze one that's not on OpenDataPhilly
Week #4 Agenda¶
Two parts:
- Part 1: More interactive data visualization: the HoloViz ecosystem
- Part 2: Getting started with raster data
Part 1: More interactive data viz¶
Recap: Data viz in Python¶
What have we learned so far¶
Matplotlib¶
- The classic, most flexible library
- Can handle geographic data well
- Overly verbose syntax, syntax is not declarative
Pandas¶
- Quick, built-in interface
- Not as many features as other libraries
seaborn¶
- Best for visualizing complex relationships between variables
- Improves matplotlib's syntax: more declarative
altair¶
- Easy, declarative syntax
- Lots of interactive features
- Complex visualizations with minimal amounts of code
We'll learn one more today...¶
A set of coordinated visualization libraries in Python¶
The motivation behind HoloViz mirrors the goals of this course¶
Proper data visualization is crucial throughout all of the steps of the data science pipeline: data wrangling, modeling, and storytelling
A significant pro:¶
GeoViews builds on HoloViews to add support for geographic data
The major con:¶
- All are relatively new
- Bokeh is the most well-tested
- HoloViews, GeoViews, hvPlot are being actively developed but are very promising
How does hvPlot fit in?¶
The hvPlot package¶
It's relatively new: officially released in February 2019
In [1]:
%%html
<center>
<blockquote class="twitter-tweet" data-lang="en"><p lang="en" dir="ltr">We are very pleased officially announce the release of hvPlot! It provides a high-level plotting API for the PyData ecosystem including <a href="https://twitter.com/pandas_dev?ref_src=twsrc%5Etfw">@pandas_dev</a>, <a href="https://twitter.com/xarray_dev?ref_src=twsrc%5Etfw">@xarray_dev</a>, <a href="https://twitter.com/dask_dev?ref_src=twsrc%5Etfw">@dask_dev</a>, <a href="https://twitter.com/geopandas?ref_src=twsrc%5Etfw">@geopandas</a> and more, generating interactive <a href="https://twitter.com/datashader?ref_src=twsrc%5Etfw">@datashader</a> and <a href="https://twitter.com/BokehPlots?ref_src=twsrc%5Etfw">@BokehPlots</a>. <a href="https://t.co/Loc5XElJUL">https://t.co/Loc5XElJUL</a></p>— HoloViews (@HoloViews) <a href="https://twitter.com/HoloViews/status/1092409050283819010?ref_src=twsrc%5Etfw">February 4, 2019</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>
</center>
Main use¶
- Quickly generate interactive plots from your data
- Seamlessly handles pandas and geopandas data
- Relies on Holoviews and Geoviews under the hood
An interface just like the pandas plot() function, but much more useful.
Let's load some data to try out¶
In [2]:
# Our usual imports
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
In [3]:
# let's load the measles data from week 2
url = "https://raw.githubusercontent.com/MUSA-550-Fall-2021/week-2/master/data/measles_incidence.csv"
measles_data_raw = pd.read_csv(url, skiprows=2, na_values='-')
In [4]:
measles_data_raw.head()
Out[4]:
Convert from wide to long formats...¶
In [5]:
measles_data = measles_data_raw.melt(id_vars=["YEAR", "WEEK"],
value_name="incidence",
var_name="state")
In [6]:
measles_data.head()
Out[6]:
Reminder: plotting with pandas¶
The default .plot() doesn't know which variables to plot.
In [7]:
fig, ax = plt.subplots(figsize=(10, 6))
measles_data.plot(ax=ax)
Out[7]:
But we can group by the year, and plot the national average each year
In [8]:
by_year = measles_data.groupby("YEAR")['incidence'].sum()
by_year.head()
Out[8]:
In [9]:
fig, ax = plt.subplots(figsize=(10, 6))
# Plot the annual average by year
by_year.plot(ax=ax)
# Add the vaccine year and label
ax.axvline(x=1963, c='k', linewidth=2)
ax.text(1963, 27000, " Vaccine introduced", ha='left', fontsize=18);
Adding interactivity with hvplot¶
Use the .hvplot() to create interactive plots.
In [10]:
# This will add the .hvplot() function to your DataFrame!
import hvplot.pandas
# Import holoviews too
import holoviews as hv
# Load bokeh
hv.extension('bokeh')
In [12]:
img = by_year.hvplot(kind='line')
img
Out[12]:
In this case, .hvplot() creates a Holoviews Curve object.
Not unlike altair Chart objects, it's an object that knows how to translate from your DataFrame data to a visualization.
In [13]:
print(img)
Many different chart types are available...¶
In [14]:
by_year.hvplot(kind='scatter')
Out[14]:
In [14]:
by_year.hvplot(kind='bar', rot=90, width=1000)
Out[14]:
Just like in altair, we can also layer chart elements together¶
Use the * operator to layer together chart elements.
Note: the same thing can be accomplished in altair, but with the + operator.
In [15]:
# The line chart of incidence vs year
incidence = by_year.hvplot(kind='line')
# Vertical line + label for vaccine year
vline = hv.VLine(1963).opts(color='black')
label = hv.Text(1963, 27000, " Vaccine introduced", halign='left')
final_chart = incidence * vline * label
final_chart
Out[15]:
We can group charts by a specific column, with automatic widget selectors¶
This is some powerful magic.
Let's calculate the annual measles incidence for each year and state:
In [16]:
by_state = measles_data.groupby(['YEAR', 'state'])['incidence'].sum()
by_state.head()
Out[16]:
Now, tell hvplot to plot produce charts for each state:
In [17]:
by_state_chart = by_state.hvplot(x="YEAR",
y="incidence",
groupby="state",
width=400,
kind="line")
by_state_chart
Out[17]:
We can select out individual charts from the set of grouped objects¶
In [20]:
PA = by_state_chart['PENNSYLVANIA'].relabel('PA')
NJ = by_state_chart['NEW JERSEY'].relabel('NJ')
Combine charts as subplots with the + operator¶
In [21]:
combined = PA + NJ
combined
Out[21]:
In [22]:
print(combined)
The charts are side-by-side by default. You can also specify the number of rows/columns explicitly.
In [25]:
# one column
combined.cols(1)
Out[25]:
We can also show overlay lines on the same plot¶
Using the by keyword:
In [26]:
states = ['NEW YORK', 'NEW JERSEY', 'CALIFORNIA', 'PENNSYLVANIA']
sub_states = by_state.loc[:, states]
In [27]:
sub_state_chart = sub_states.hvplot(x='YEAR',
y='incidence',
by='state',
kind='line')
sub_state_chart * vline
Out[27]:
We can also show faceted plots¶
Just like in altair, when we used the alt.Chart().facet(column='state') syntax
Below, we specify the state column should be mapped to each column:
In [29]:
img = sub_states.hvplot(x="YEAR",
y='incidence',
col="state",
kind="line",
rot=90,
frame_width=200)
img * vline
Out[29]:
Functions for each kind of chart type are available too¶
In [25]:
by_state.hvplot
For example, we could plot a bar chart for these four states¶
In [27]:
by_state.loc[1960:1970, states].hvplot.bar(x='YEAR',
y='incidence',
by='state', rot=90)
Out[27]:
Change bar() to line() and we get the same thing as before.
In [30]:
by_state.loc[1960:1970, states].hvplot.line(x='YEAR',
y='incidence',
by='state', rot=90)
Out[30]:
Customizing charts¶
See the help message for explicit hvplot functions:
In [32]:
by_state.hvplot?
In [34]:
by_state.hvplot.line?
Heatmaps are available too...¶
Can we reproduce the WSJ measles heatmap that we made in altair in week 2?
Use the help function:
In [35]:
measles_data.hvplot.heatmap?
Two methods:¶
We want to plot 'YEAR' on the x axis, 'state' on the y axis, and specify 'incidence' as the values begin plotted in each heatmap bin.
- You can use the
by_statedata frame which has already summed over weeks for each state - Use the original, tidy data (
measles_data) with columns for state, week, year, and incidence- you will need to use the
reduce_functionkeyword to sum over weeks
- you will need to use the
In [36]:
by_state
Out[36]:
In [38]:
# METHOD #1: just plot the incidence
heatmap = by_state.hvplot.heatmap(
x="YEAR",
y="state",
C="incidence",
cmap="viridis",
height=500,
width=1000,
flip_yaxis=True,
rot=90,
)
heatmap.redim(
state="State", YEAR="Year",
)
Out[38]:
In [ ]:
In [40]:
## METHOD 2: hvplot does the aggregation
heatmap = measles_data.hvplot.heatmap(
x="YEAR",
y="state",
C="incidence",
cmap='viridis',
reduce_function=np.sum,
height=500,
width=1000,
flip_yaxis=True,
rot=90,
)
heatmap.redim(state="State", YEAR="Year")
Out[40]:
Just like altair: save the file as html¶
In [41]:
import hvplot
hvplot.save(heatmap, 'measles.html')
In [42]:
# load the html file and display it
from IPython.display import HTML
HTML('measles.html')
Out[42]:
Let's load the penguins data set from week 2
In [43]:
url = "https://raw.githubusercontent.com/MUSA-550-Fall-2021/week-2/master/data/penguins.csv"
penguins = pd.read_csv(url)
In [44]:
penguins.head()
Out[44]:
Use the hvplot.scatter_matrix() function:
In [35]:
penguins.hvplot.scatter?
In [45]:
columns = ['flipper_length_mm',
'bill_length_mm',
'body_mass_g',
'species']
hvplot.scatter_matrix(penguins[columns], c='species')
Out[45]:
Note the "box select" and "lasso" features on the tool bar for interactions
Recap: altair vs hvplot¶
- Both use a declarative syntax (altair more so than hvplot)
- Users of ggplot might be more familiar with altair's syntax
- hvplot integrates directly into pandas dataframes via the
.hvplot()function
- Both have support for cross-filtering and interactions
- Both can be incorporated into web-based dashboard via HTML (later in course)
- hvplot has better support for large data (later in course)
It's largely up to you which one you feel is easier to use¶
hvplot can also be used with geopandas!¶
Let's load some geographic data for countries:
In [46]:
import geopandas as gpd
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
In [47]:
world.head()
Out[47]:
Plotting with just geopandas¶
In [48]:
fig, ax = plt.subplots(figsize=(10,10))
world.plot(column='gdp_md_est', ax=ax)
ax.set_axis_off()
Now with hvplot¶
In [49]:
world.hvplot.polygons?
In [51]:
# Can also just do world.hvplot()
world.hvplot.polygons(c='gdp_md_est',
geo=True,
frame_height=400)
Out[51]:
Let's try it on our median assessment values per neighborhood from last week¶
In [45]:
import geopandas as gpd
In [52]:
# Load the data
url = "https://raw.githubusercontent.com/MUSA-550-Fall-2021/week-3/master/data/opa_residential.csv"
data = pd.read_csv(url)
# Create the Point() objects
data['Coordinates'] = gpd.points_from_xy(data['lng'], data['lat'])
# Create the GeoDataFrame
data = gpd.GeoDataFrame(data, geometry='Coordinates', crs="EPSG:4326")
In [53]:
# load the Zillow data from GitHub
url = "https://raw.githubusercontent.com/MUSA-550-Fall-2020/week-3/master/data/zillow_neighborhoods.geojson"
zillow = gpd.read_file(url)
In [54]:
# Important: Make sure the CRS match
data = data.to_crs(zillow.crs)
# perform the spatial join
data = gpd.sjoin(data, zillow, op='within', how='left')
In [55]:
# Calculate the median market value per Zillow neighborhood
median_values = data.groupby('ZillowName', as_index=False)['market_value'].median()
# Merge median values with the Zillow geometries
median_values = zillow.merge(median_values, on='ZillowName')
print(type(median_values))
In [56]:
median_values.head()
Out[56]:
In [57]:
median_values.crs
Out[57]:
In [52]:
# pass arguments directly to hvplot()
# and it recognizes polygons automatically
median_values.hvplot(c='market_value',
frame_width=600,
frame_height=500,
geo=True,
cmap='viridis',
hover_cols=['ZillowName'])
Out[52]:
Important: geo=True assumes EPSG:4326¶
If you specify geo=True, the data needs to be in typical lat/lng CRS. If not, you can use the crs keyword to specify the type of CRS your data is in.
In [59]:
median_values_3857 = median_values.to_crs(epsg=3857)
In [60]:
median_values_3857.crs
Out[60]:
In [61]:
median_values_3857.hvplot(c='market_value',
frame_width=600,
frame_height=500,
geo=True,
crs=3857, # NEW: specify the CRS
cmap='viridis',
hover_cols=['ZillowName'])
Out[61]:
Now we can take advantage of GeoViews¶
Let's add a tile source underneath the choropleth map
In [62]:
import geoviews as gv
import geoviews.tile_sources as gvts
In [57]:
%%opts WMTS [width=800, height=800, xaxis=None, yaxis=None]
choro = median_values.hvplot(c='market_value',
width=500,
height=400,
alpha=0.5,
geo=True,
cmap='viridis',
hover_cols=['ZillowName'])
gvts.ESRI * choro
Out[57]:
Many of the most common tile sources are available..¶
In [64]:
print(type(gvts.ESRI))
In [65]:
%%opts WMTS [width=200, height=200, xaxis=None, yaxis=None]
(gvts.OSM +
gvts.Wikipedia +
gvts.StamenToner +
gvts.EsriNatGeo +
gvts.EsriImagery +
gvts.EsriUSATopo +
gvts.EsriTerrain +
gvts.CartoDark).cols(4)
Out[65]:
Note: we've used the %%opts cell magic to apply syling options to any charts generated in the cell.
See the documentation guide on customizations for more details.
What about interactive hex bins?¶
You can do it with hvplot! Sort of.
Step 1: Extract out the x/y values of the data¶
- Let's add them as new columns into the data frame
- Remember, you can the use "x" and "y" attributes of the "geometry" column.
In [66]:
data['x'] = data.geometry.x
data['y'] = data.geometry.y
In [61]:
data.head()
Out[61]:
Step 2: Create a new DataFrame with only the columns we need¶
- In this case, we'll use the
xandycoordinate columns and the associatedmarket_valuecolumn - Important: hvplot's hexbin won't work if the input data is a GeoDataFrame, e.g., still has a
geometrycolumn- By doing the column selection, we convert the data into a regular pandas DataFrame
In [68]:
subdata = data[['x', 'y', 'market_value']]
In [69]:
type(subdata)
Out[69]:
Step 3: Plot with the hexbin function¶
- Similar syntax to matplotlib's hexbin() function
- Specify:
- The x/y coordinates,
- An optional
Ccolumn to aggregate for each bin (raw counts are shown if not provided) - A
reduce_functionthat determines how to aggregate theCcolumn
In [64]:
data.hvplot.hexbin?
In [65]:
subdata.head()
Out[65]:
In [66]:
subdata.hvplot.hexbin(x='x',
y='y',
C='market_value',
reduce_function=np.median,
logz=True,
geo=True,
gridsize=40,
cmap='viridis')
Out[66]:
Not the prettiest but it gets the job done for some quick exploratory analysis!
Documentation references¶
- Hvplot user guide
- HoloViz tutorial: introduction to the HoloViz ecosystem
- HoloViews user guide and gallery
- GeoViews user guide and gallery
Some very cool examples available in the galleries
That's it!¶
- Homework #2 due on Wednesday
- Office hours tomorrow evening, 6-8pm (sign up on Canvas)
- Wednesday we'll dive into raster data analysis
In [ ]: