Skip to article frontmatterSkip to article content

In this lesson, we’ll learn about geospatial operations such as dissolves, intersections, and spatial joins. We’ll also relate tabular data to geospatial data and learn new ways of combining different types of datasets. By the end of this lesson, students will be able to:

  • Apply the dissolve method to aggregate GeoDataFrame columns.

  • Apply merge to join tabular data and sjoin to join geospatial data.

  • Apply the clip methods to perform bounding box queries.

import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd

We’ll continue using the same world countries dataset today.

columns = ["POP_EST", "GDP_MD", "CONTINENT", "SUBREGION", "geometry"]
countries = gpd.read_file("ne_110m_admin_0_countries.shp").set_index("NAME")[columns]
countries

Dissolve: groupby for geospatial data

How can we plot the total population for each continent (rather than per country)? In pandas, to perform a computation for each group, we would typically use a groupby!

populations = countries.groupby("CONTINENT")["POP_EST"].sum()
populations.plot()

Clearly, this is not a map. For geospatial data, we need to use a GeoDataFrame method called dissolve. Dissolve not only takes care of aggregation, but it also performs a geospatial union to combine all the countries within a single continent into a larger geometric shape.

Combining tabular and geospatial data

When we first learned about data frames, we introduced a dataset about earthquakes around the world. How might we plot these earthquakes?

earthquakes = gpd.read_file("earthquakes.csv").set_index("id")
earthquakes

Unfortunately, even though the earthquakes dataset has latitude and longitude columns, it is not strictly in a geospatial data format like a shapefile. It’s better to read this dataset using plain pandas and then work on adding a geometry column from the latitude-longitude pairs.

earthquakes = pd.read_csv("earthquakes.csv").set_index("id")
earthquakes = gpd.GeoDataFrame(
    earthquakes,
    # crs="EPSG:4326" specifies WGS84 or GPS coordinate system, see https://epsg.io/4326
    geometry=gpd.points_from_xy(earthquakes["longitude"], earthquakes["latitude"], crs="EPSG:4326")
)
earthquakes

Now, we can plot the earthquakes on a map!

fig, ax = plt.subplots(figsize=(13, 5))
countries.plot(ax=ax, color="#EEE")
earthquakes.plot(ax=ax, column="magnitude", markersize=0.1, legend=True)
ax.set(title="Earthquakes between July 27, 2016 and August 25, 2016")
ax.set_axis_off()

But what if we want to only plot the earthquakes that occurred over land in North America? The following map shows the earthquakes in Washington atop a background of North America, but it’s pretty tedious trying to specify every single place name one-by-one.

fig, ax = plt.subplots()
na_countries = countries[countries["CONTINENT"] == "North America"]
na_countries.plot(ax=ax, color="#EEE")
earthquakes[earthquakes["name"] == "Washington"].plot(ax=ax, column="magnitude", markersize=0.1)
ax.set(title="Earthquakes in Washington")
ax.set_axis_off()

Bounding boxes

A bounding box is one way to specify an intersection query. Say we wanted to show all the earthquakes over the map area covered by North America. We can see the bounding boxes for every country in North America through the bounds field.

na_countries.bounds

The total_bounds field provides the bounding box for the entire dataset.

na_countries.total_bounds

We can use the total bounds as input to the clip method to keep only the rows where the geometry falls within the total bounds.

earthquakes_na_bounds = earthquakes.clip(na_countries.total_bounds)
earthquakes_na_bounds
fig, ax = plt.subplots()
na_countries = countries[countries["CONTINENT"] == "North America"]
na_countries.plot(ax=ax, color="#EEE")
earthquakes_na_bounds.plot(ax=ax, column="magnitude", markersize=0.1)
ax.set(title="Earthquakes in North America (over land or sea)")
ax.set_axis_off()

Spatial join

But what if we only want to plot the earthquakes that appeared over land—not over sea? A spatial join allows us to specify geospatial operations to link two datasets so that we can find all the earthquakes that were located in a North American country. This allows us to provide more precise intersection queries.

earthquakes_na_countries = earthquakes.sjoin(na_countries)
earthquakes_na_countries

The columns to the left of geometry are from the earthquakes dataset, while the columns to the right of geometry are linked-up from the na_countries dataset. Complete this sentence to describe how sjoin computed the above result.

For every row in the earthquakes dataset, sjoin _________________________

By default, sjoin uses the keyword argument how="inner". We can also customize the join result with how="left" or how="right" too. Can you explain how the results differ?

earthquakes.sjoin(na_countries, how="inner")

Finally, let’s plot all the earthquakes that occurred over land in North America.

fig, ax = plt.subplots()
na_countries = countries[countries["CONTINENT"] == "North America"]
na_countries.plot(ax=ax, color="#EEE")
earthquakes_na_countries.plot(ax=ax, column="magnitude", markersize=0.1)
ax.set(title="Earthquakes in North America (over land)")
ax.set_axis_off()

Attribute joins

The idea of joining two datasets isn’t exclusive to geospatial data. In fact, in many data-centric contexts we may need to combine multiple datasets. Just as sjoin combines two geospatial datasets on their geometry columns using geometric intersection, merge combines two tabular datasets on specific columns using exact == matches.

Consider the following two datasets of movies and directors. Not all directors have a movie listed, and not all movies have a corresponding director in the dataset.

movies = pd.DataFrame([
    {"movie_id": 51, "movie_name": "Lady Bird", "year": 2017, "director_id": 23},
    {"movie_id": 47, "movie_name": "Grand Budapest Hotel", "year": 2014, "director_id": 16},
    {"movie_id": 103, "movie_name": "Parasite", "year": 2019, "director_id": 14},
    {"movie_id": 34, "movie_name": "Frozen", "year": 2013, "director_id": 18},
    {"movie_id": 37, "movie_name": "Moonrise Kingdom", "year": 2012, "director_id": 16},
])
movies
directors = pd.DataFrame([
    {"director_id": 14, "director_name": "Bong Joon Ho"},
    {"director_id": 23, "director_name": "Greta Gerwig"},
    {"director_id": 16, "director_name": "Wes Anderson"},
    {"director_id": 21, "director_name": "Quentin Tarantino"},
    {"director_id": 27, "director_name": "Kathryn Bigelow"},
])
directors

Before running the following cell, predict the output if the default attribute join type is how="inner". Visualize the procedure using PandasTutor. What would happen if we changed the join type?

movies.merge(directors, on="director_id")

For these two datasets, the column names are the same so we can use the on keyword argument. If the column names are not the same, specify the column name in the left with left_on and the column name in the right with right_on.

Interactive maps

Geopandas supports interactive maps with the folium library.

!pip install -q folium mapclassify

Generating an interactive map is as simple as generating a static map: instead of calling plot, call explore.

earthquakes.explore(column="magnitude")
Lectures
Geospatial Data Frames
Lectures
Model Evaluation