Week 3: Geospatial Data Analysis and GeoPandas¶
Sep 15, 2020
Housekeeping¶
- Homework #2 due a week of Thursday (9/24)
- 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
Agenda for Week #3¶
- Vector data and introduction to GeoPandas
- Spatial relationships and joins
- Visualization for geospatial data
- Demo: 311 requests by neighborhood in Philadelphia
- Exercise: Property assessments by neighborhood
Now, on to geospatial analysis...¶
# Let's setup the imports we'll need first
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
%matplotlib inline
Vector Data¶
- Vector refers to discrete geometric entities
- The Open Geospatial Consortium has standardized a set of simple features
- Includes points, lines, and polygons
A couple of terminology notes¶
- A feature refers to both the geometry and attributes of specific piece of vector data
- A feature collection is a list, or collection, of features
Both terms are very common in Python geospatial software.
Common formats for vector datasets¶
A shapefile¶
Actually several files with the same common prefix
Mandatory files:
- .shp: the file containing the geometries
- .shx: the file that indexes the geometry
- .dbf: tabular data format storing the attributes for each geometry
And many optional files for documentation, projection information, etc.
Let's take a look at an example shapefile:
We'll use the ls command to list out all of the files in an example shapefile in the data/ folder
%ls "data/ne_110m_admin_0_countries/"
The GeoJSON file¶
- Stores simple features in a JSON format
- Arose due to the prevalence of the JSON format, especially on the web
Additional GeoJSON resources and tools:¶
- GitHub lets you view GeoJSON files natively
- http://geojson.io provides interactive creation and viewing of small GeoJSON files
GitHub example from the data/ directory: Philadelphia ZIP Codes
Working with vector data in Python: GeoPandas¶
geopandas provides a simple, intuitive for the main types of geospatial vector file formats
Example: Let's load a shape file of countries in the world...
import geopandas as gpd
We can use the read_file() function to read shapefiles and GeoJSON files.
# Read the shape file, giving the name of the directory
countries = gpd.read_file("./data/ne_110m_admin_0_countries")
countries.head()
type(countries)
What's a GeoDataFrame?¶
Just like a DataFrame but with a geometry column
# Print out the first 10 entires of the "geometry" column
countries.geometry.head(n=10)
We can still leverage the power of pandas...¶
Calculate the total world population:
countries['pop_est'].sum()/1e9 # In billions
Calculate the total population on each continent:
grouped = countries.groupby('continent')
grouped
Remember: the
groupby() does not return a DataFrame — you need to call sum(), mean(), etc, or apply() a function.
# Sum population on each continent
pop_by_continent = grouped['pop_est'].sum()
# Sort values
pop_by_continent.sort_values(ascending=False, inplace=True)
# Output sorted values from cell
pop_by_continent
Filter the data frame based on a boolean selection:
# Is the country name USA?
is_USA = countries['name']=='United States of America'
# Get the row with USA
USA = countries.loc[is_USA]
USA
An aside: the squeeze() function¶
It does just one it sounds like: if you have a DataFrame with only one row, it will "squeeze" the row dimension by removing it, returning just a Series object:
# Squeeze
USA = USA.squeeze()
# Print out the type
print("The type of USA is: ", type(USA))
# Output
USA
The simple features (Lines, Points, Polygons) are implemented by the shapely library
type(USA.geometry)
Jupyter notebook renders shapely geometries automatically:
# a mini USA
USA.geometry
How does geopandas handle coordinate systems and map projections?¶
Coordinate Reference Systems¶
A coordinate reference system (CRS) relates the position of a geometry object on the spherical earth to its two-dimensional coordinates.
A GeoDataFrame or GeoSeries has a .crs attribute which specifies the coordinate reference system.
countries.crs
- EPSG 4326 is known as WGS 84 where
xandyare longitude and latitude. - It is is the default coordinate system for GPS systems.
- It's also known as Plate Carrée
How to plot all of the geometries at once?¶
Use the plot() function to get a quick and dirty plot of all of the geometry features.
Note: the plot() returns the current maplotlib axes, allowing you to format the chart after plotting.
# Create a figure and axes
fig, ax = plt.subplots(figsize=(10, 6))
# Plot the countries on our axes
ax = countries.plot(ax=ax, facecolor="none", edgecolor="black")
# Add a title
ax.set_title("Equirectangular Projection")
What's going on under the hood?¶
matplotlib and cartopy are combined to make geo-aware plots
import cartopy.crs as ccrs
# Initialize the EPSG:4326 CRS object
wgs84 = ccrs.PlateCarree()
# Create a geo-aware axes using the "projetion" keyword
fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={"projection": wgs84})
# Print out the type of the axes for info purposes
print("The type of the axes is: ", type(ax))
# Add the geometry shapes
ax.add_geometries(countries["geometry"], crs=wgs84, facecolor="none", edgecolor="black")
# Add a title
ax.set_title("Equirectangular Projection");
See the Geopandas documentation for more examples: Plotting with CartoPy and GeoPandas
Can we convert to other geometries?¶
Use the to_crs() function. The most well-known projections can be specified by their EPSG code.
Geopandas documentation on re-projecting: Managing Projections
Let's convert to the Mercator projection¶
# Remove Antartica, as the Mercator projection
# cannot deal with the poles
no_antartica = countries[(countries['name'] != "Antarctica")]
# Two ways to specify the EPSG code
countries_mercator = no_antartica.to_crs(epsg=3395)
# Alternatively:
# countries_mercator = no_antartica.to_crs("EPSG:3395")
countries_mercator.head()
Note: the magnitude of the values in the geometry column changed! A quick and easy way to tell if the re-projection worked properly!
Now let's plot it¶
The easy way...with geopandas built-in plot() function
# Initialize the figure and axes
fig, ax = plt.subplots(figsize=(10, 6))
# Use built-in plot() of the GeoDataFrame
ax = countries_mercator.plot(ax=ax, facecolor="none",
edgecolor="black")
# Add a title
ax.set_title("Mercator Projection");
The harder way...¶
With cartopy and matplotlib
# Initialize the CRS object
crs = ccrs.epsg(3395) # or crs = ccrs.Mercator()
# Create a geo-aware axes using the "projetion" keyword
fig, ax = plt.subplots(figsize=(10, 6),
subplot_kw={"projection": crs})
# Add the geometry shapes
ax.add_geometries(countries_mercator['geometry'],
crs=crs, facecolor='none', edgecolor='black')
# Add a title
ax.set_title('Mercator Projection');
When to use the 'harder' way?¶
When you need more customizable, advanced plots.
Nearly anything that matplotlib can do can be plotted on a cartopy GeoAxes. Plotting directly with matplotlib allows you to take full advantage of matplotlib's functionality.
So which CRS is best?¶
- For city-based data, usually Web Mercator (EPSG=3857) is best
- Can also use a CRS specific to individual states, e.g., PA State Plane EPSG=2272
Let's load the city limits for Philadelphia¶
We'll use the provided City_Limits shape file in the data/ folder
city_limits = gpd.read_file('./data/City_Limits')
city_limits
What's the CRS?¶
Use the .crs attribute to find out!
city_limits.crs
# Create our figure and axes
fig, ax = plt.subplots(figsize=(5, 5))
# Plot
city_limits.plot(ax=ax, facecolor="none", edgecolor="black")
# Format
ax.set_title("Equirectangular")
ax.set_axis_off() # This will remove the axes completely
ax.set_aspect("equal") # This forces an equal aspect ratio
This is not what Philadelphia looks like..
Let's try EPSG=3857 instead:
# Create the figure
fig, ax = plt.subplots(figsize=(5, 5))
# Convert to EPSG:3857
city_limits_3857 = city_limits.to_crs(epsg=3857)
# Plot and format
city_limits_3857.plot(ax=ax, facecolor="none", edgecolor="black")
ax.set_title("Web Mercator")
ax.set_axis_off()
ax.set_aspect("equal");
Important: the equirectangular CRS (EPSG=4326) is often used by default and will make cities appear wider and flatter than they really are
Saving GeoDataFrames¶
Use the to_file() function and specify the driver.
# ESRI shape file
city_limits_3857.to_file("./data/city_limits_3857",
driver='ESRI Shapefile')
# GeoJSON is also an option
city_limits_3857.to_file("./data/city_limits_3857.geojson",
driver='GeoJSON')
How about as a CSV file?¶
Yes, but reading requires more work...
# save a csv file
city_limits_3857.to_csv("./data/city_limits_3857.csv", index=False)
df = pd.read_csv("./data/city_limits_3857.csv")
df.head()
Looks similar...¶
But, the "geometry" column is just stored as a string...it's not a shapely Polygon
type(df.geometry.iloc[0])
Use shapely to parse the string version of the polygons¶
from shapely import wkt
# wkt.loads will convert from string to Polygon object
df['geometry'] = df['geometry'].apply(wkt.loads)
type(df.geometry.iloc[0])
Converting from a DataFrame to a GeoDataFrame¶
We can initialize the GeoDataFrame directly from a DataFrame but we need to specify two things:
- The name of the "geometry" column
- The CRS of the "geometry" column
In this case, the geometry column was saved in Web Mercator EPSG=3857
# Make specifying the name of the geometry column and CRS
gdf= gpd.GeoDataFrame(df, geometry='geometry', crs="EPSG:3857")
# Now plot
fig, ax = plt.subplots(figsize=(5,5))
ax = gdf.plot(ax=ax, facecolor='none', edgecolor='black')
ax.set_axis_off()
ax.set_aspect("equal")
Let's convert back to 4326 and plot¶
The tilt should be a bit more obvious now...
ax = gdf.to_crs(epsg=4326).plot(facecolor='none', edgecolor='black')
ax.set_axis_off()
ax.set_aspect("equal")
Note¶
- I didn't use
plt.subplots()here to create a figure/axes – I let geopandas automatically make one - I've chained together the
to_crs()and.plot()functions in one line - The
.plot()function returns the axes object that geopandas used to plot — this lets you customizes the axes after plotting
So, when should you use GeoPandas?¶
- For exploratory data analysis and visualization, including in Jupyter notebooks
- Pre-processing data to be fed into a desktop GIS program
- For compact, readable, and reproducible code
- If you’re comfortable with Pandas and/or R data frames.
When it may not be the best tool:¶
- For polished multilayer map creation — one option is to use a desktop GIS like QGIS.
- If you need very high performance — geopandas can be slow compared to other GIS software.
Spatial Relationships and Joins¶
# Load some cities data
cities = gpd.read_file("./data/ne_110m_populated_places")
(Image by Krauss, CC BY-SA 3.0)
All of these operations are available as functions of a GeoDataFrame.
# Select the Point representing New York City
new_york = cities.loc[cities['name'] == 'New York', 'geometry'].squeeze()
new_york
type(new_york)
countries.contains(new_york)
# Find the country that contains New York
countries.loc[countries.contains(new_york)]
# Get the geometry column of the country containing NYC
USA = countries.loc[countries.contains(new_york), 'geometry'].squeeze()
USA
Note¶
The .loc[] function can take the index selector as the first argument, and the name of a column as a second argument (separated by a comma)
type(USA)
# Is New York within the USA?
new_york.within(USA)
Reference¶
The different functions for checking spatial relationships:
equalscontainscrossesdisjointintersects
overlapstoucheswithincovers
See the shapely documentation for an overview of these methods.
The spatial join¶
SPATIAL JOIN = merging attributes from two geometry layers based on their spatial relationship
Different parts of this operations:
- The GeoDataFrame to which we want add information
- The GeoDataFrame that contains the information we want to add
- The spatial relationship we want to use to match both datasets (intersects, contains, within)
- The type of join: left or inner join
In this case, we want to join the cities dataframe, containing Point geometries, with the information of the countries dataframe, containing Polygon geometries.
To match cities with countries, we'll use the within spatial relationship.
The geopandas.sjoin() function performs this operation:
joined = gpd.sjoin(cities, countries, op='within', how='left')
joined.head()
# Extract Italy
italy = countries.loc[countries['name']=='Italy']
# Plot
fig, ax = plt.subplots(figsize=(8,8))
italy.plot(ax=ax, facecolor='none', edgecolor='black')
ax.set_axis_off()
ax.set_aspect("equal")
# Plot the first city in the joined data frame (Vatican City)
# Use the same axes by passing in the ax=ax keyword
ax = joined.iloc[:1].plot(ax=ax, color='red')
Spatial overlay operation¶
We can also perform the join() operation on the geometries rather than just combining attributes.
The overlay() function combines geometries, e.g. by taking the intersection of the geometries.
africa = countries[countries['continent'] == 'Africa']
## What crs?
africa.crs
# Let's transform to a CRS that uses meters
# instead of degrees (EPSG=3857)
africa = africa.to_crs(epsg=3857)
africa.crs
fig, ax = plt.subplots(figsize=(8,8))
africa.plot(ax=ax, facecolor='#b9f2b1')
ax.set_axis_off()
ax.set_aspect("equal")
# Important CRS needs to match!
cities = cities.to_crs(epsg=3857)
# Create a copy of the GeoDataFrame
buffered_cities = cities.copy()
# Add a buffer region of 250 km around all cities
buffered_cities['geometry'] = buffered_cities.buffer(250e3)
Plot the difference of the two geometries¶
fig, ax = plt.subplots(figsize=(8, 8))
# Calculate the difference of the geometry sets
diff = gpd.overlay(africa, buffered_cities, how='difference')
# Plot
diff.plot(facecolor="#b9f2b1", ax=ax)
ax.set_axis_off()
ax.set_aspect("equal")
# Data attributes are the same as the first data frame (africa)
# with an updated geometry column
diff.head()
Plot the intersection of the two geometries¶
fig, ax = plt.subplots(figsize=(8, 8))
# The intersection of the geometry sets
intersection = gpd.overlay(africa, buffered_cities, how='intersection')
# Plot
intersection.plot(ax=ax, facecolor="#b9f2b1")
ax.set_axis_off()
ax.set_aspect("equal")
Recap:¶
- Spatial join: merge attributes from one data frame to another based on the spatial relationship
- Spatial overlay: creating new geometries based on spatial operation between both data frames (and not combining attributes of both data frames)
Putting it all together: 311 requests in 2020¶
Load 311 requests in Philadelphia from the data/ directory.
Source: OpenDataPhilly
# Load the data from a CSV file into a pandas DataFrame
requests = pd.read_csv('./data/public_cases_fc_2020.csv')
print("number of requests = ", len(requests))
requests.head()
First, convert to a GeoDataFrame¶
Remove the requests missing lat/lon coordinates
requests = requests.dropna(subset=['lat', 'lon'])
Create Point objects for each lat and lon combination.
We can use the helper utility function: geopandas.points_from_xy()
requests['Coordinates'] = gpd.points_from_xy(requests['lon'], requests['lat'])
requests['Coordinates'].head()
Now, convert to a GeoDataFrame.
Important
- Don't forget to set the CRS manually!
- The CRS you specify when creating a GeoDataFrame should tell geopandas what the coordinate system the input data is in.
- Usually you will be reading lat/lng coordinates, and will need to specify the crs as EPSG code 4326
- You should specify the crs as a string using the syntax:
ESPG:4326
Since we're only using a few EPSG codes in this course, you can usually tell what the CRS is by looking at the values in the Point() objects.
Philadelphia has a latitude of about 40 deg and longitude of about -75 deg.
Our data must be in the usual lat/lng EPSG=4326.
requests = gpd.GeoDataFrame(requests, geometry="Coordinates", crs="EPSG:4326")
Next, identify the top 20 most common requests¶
Group by the service name and calculate the size of each group:
service_types = requests.groupby('service_name').size()
Sort by the number (in descending order):
service_types = service_types.sort_values(ascending=False)
Slice the data to take the first 20 elements:
top20 = service_types.iloc[:20]
top20
Let's trim to only the trash-related requests¶
trash_requests = requests.loc[
requests["service_name"] == "Rubbish/Recyclable Material Collection"
]
print("The nuumber of trash-related requests = ", len(trash_requests))
Trash collection has been a BIG issue in Philadelphia recently¶
See for example, this article in the Philadelphia Inquirer
Let's plot the monthly totals for 2020¶
# Convert the requested datetime to a column of Datetime objects
trash_requests['requested_datetime'] = pd.to_datetime(trash_requests['requested_datetime'])
# Use the .dt attribute to extract out the month name
trash_requests['month'] = trash_requests['requested_datetime'].dt.month_name()
totals_by_week = trash_requests.groupby("month", as_index=False).size()
totals_by_week.head()
Note: I've used the as_index=False syntax here¶
This will force the size() function to return a DataFrame instead of having the month column as the index of the resulted groupby operation.
It saves us from having to do the .reset_index() function call after running the .size() function.
Plot a bar chart with seaborn¶
For making static bar charts with Python, seaborn's sns.barplot() is the best option
import seaborn as sns
# Initialize figure/axes
fig, ax = plt.subplots(figsize=(12, 6))
# Plot!
sns.barplot(
x="month",
y="size",
data=totals_by_week,
color="#2176d2",
ax=ax,
order=["January", "February", "March", "April", "May", "June", "July", "August", "September"],
);
Example: Improving the aesthetics of matplotlib¶
The trend is clear in the previous chart, but can we do a better job with the aesthetics? Yes!
For reference, here is a common way to clean up charts in matplotlib:
# Initialize figure/axes
fig, ax = plt.subplots(figsize=(12, 6))
# Plot!
sns.barplot(
x="month",
y="size",
data=totals_by_week,
color="#2176d2",
ax=ax,
order=["January", "February", "March", "April", "May", "June", "July", "August", "September"],
zorder=999 # Make sure the bar charts are on top of the grid
)
# Remove x/y axis labels
ax.set_xlabel("")
ax.set_ylabel("")
# Format the ytick labels to use a comma and no decimal places
ax.set_yticklabels([f"{yval:,.0f}" for yval in ax.get_yticks()] )
# Add a grid backgrou d
ax.grid(True, axis='y')
# Remove the top and right axes lines
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# Add a title
ax.set_title("Philadelphia's Trash-Related 311 Requests in 2020", weight='bold', fontsize=16);
