Adding Basemaps In Python With Contextily

If you work with geospatial data in Python, you most likely are familiar with the fantastic GeoPandas library. GeoPandas leverages the power of Maplotlib to enable users to make maps of their data. However, until recently, it has not been easy to add basemaps to these maps. Basemaps are the contextual map data, like Google Maps, on top of which geospatial data are often displayed.

The new Python library contextily, which stands for contextual map tiles, now makes it possible and relatively straightforward to add map tile basemaps to Geopandas maps or to save map tiles as geospatial raster files. Below we walk through a few common workflows for doing this.

We start by importing the libraries we will use.  

%matplotlib inline

import pandas as pd

import geopandas as gpd

import contextily as cx

import matplotlib.pyplot as plt

This assumes we have already installed in our Python environment the following libraries

• pandas
• geopandas (and all dependancies)
• matplotlib
• contextily
• descartes

Installing geopandas can be tricky - so please read the installation instructions on the GeoPandas website!

READ DATA INTO A GEOPANDAS GEODATAFRAME

Let's fetch the census places data to map. Census places include cities and other populated places. Here we fetch the 2019 cartographic boundary (cb_) file of California (06) places. We read the data in from this URL to a GeoPandas GeoDataFrame named places.

ca_places_url = "https://www2.census.gov/geo/tiger/GENZ2019/shp/cb_2019_06_place_500k.zip"

places = gpd.read_file(ca_places_url)

Now, we can use the geodataframe plot method to make a quick map.

places.plot();

We can see the outlines of CA cities in the plotted map. Now, let's take a look at the data. It looks just like a pandas dataframe, and it is. So all pandas dataframe methods can be used with this GeoDataFrame. What is special is that GeoPandas adds methods and attributes specific to geospatial vector data. These apply to the geometry column.

places.head()

Like a pandas datafrmae, we can subset the geodataframe by selecting a row or rows by place name. Let's select the city of Berkeley, CA and map it.

berkeley = places[places['NAME']=='Berkeley']

berkeley.plot();

USE CONTEXTILY TO ADD A BASEMAP

Above we can see the map of the boundary of the city of Berkeley, CA. The axis labels display the longitude and latitude coordinates for the bounding extent of the city.

Let's use contextily in it's most simple form to add a basemap to provide the geographic context for Berkeley.

ax = berkeley.to_crs('EPSG:3857').plot(figsize=(9, 9))

cx.add_basemap(ax)

There are a few important things to note about the above code.

- We use matplotlib to define the plot canvas as ax.

- We then add the contextily basemap to the map with the code cx.add_basemap(ax)

Additionally, we dynamically transform the coordinate reference system, or CRS, of the Berkeley geodataframe from geographic lat/lon coordinates to web mercator using the method to_crs('EPSG:3857')Web mercator is the default CRS used by all web map tilesets. It is referenced by the code EPSG:3857 where EPSG stands for the initials of the organization that created these codes (the European Petroleum Survey Group).

Let's clean up the map by adding some code to change the symbology of the Berkeley city boundary. This will highlight the value of adding a basemap.

First, let's map the boundary without a fill color.

berkeley.plot(edgecolor="red", facecolor="none");

Let's build on those symbology options and add the contextily basemap.

ax = berkeley.to_crs('EPSG:3857').plot(edgecolor="red", 

                                       facecolor="none",  

                                       linewidth=2,

                                       figsize=(9, 9)

                                      )

cx.add_basemap(ax)

MAPPING POINT DATA

Let's expand on this example by mapping a point dataset of BART station locations.

First we fetch these data from a D-Lab web mapping tutorial.

bart_url = 'https://raw.githubusercontent.com/dlab-berkeley/Geospatial-Fundamentals-in-R-with-sf/master/data/bart.csv'

bart = pd.read_csv(bart_url)

Take a look at the data:

bart.head()

CONVERTING POINT DATA IN A DATAFRAME TO GEOSPATIAL DATA IN A GEODATAFRAME

Because these data are in a CSV file we read them into a Pandas DataFrame.

In order to map these data we need to convert these data to a GeoPandas GeoDataFame. To do this, we need to specify:

- the data, here the geodataframe bart

- the coordinate data, here the columns bart['X'] and bart['Y']

- the CRS of the bart coordinate data, here EPSG:4326

The CRS code 'EPSG:4326' stands for the World Geodetic System of 1984, or WGS84. This is the most commonly used CRS for geographic (lat/lon) coordinate data.

#Convert the DataFrame to a GeoDataFrame. 

bart_gdf = gpd.GeoDataFrame(bart, geometry=gpd.points_from_xy(bart['X'], bart['Y']), crs='EPSG:4326') 

# and take a look

bart_gdf.plot();

Now that we have the BART data in a geodataframe we can use the same commands as we did above to map it with a contextily basemap.

ax = bart_gdf.to_crs('EPSG:3857').plot(figsize=(9, 9))

cx.add_basemap(ax)

We have the full range of matplotlib style options to enhance the map, a few of which are shown in the example below.

ax = bart_gdf.to_crs('EPSG:3857').plot(color="red",

                                    edgecolor="black",

                                    markersize=50, 

                                    figsize=(9, 9))

ax.set_title('Bay Area Bart Stations')

cx.add_basemap(ax)

CHANGING THE BASEMAP

By default contextiley returns maptiles from the OpenStreetmap Mapnik basemap. However, there are other available tilesets from different providers. These tilesets are stored in the contextily cx.providers dictionary.

That's a large dictionary and you can view it. Alternatively, and more simply, you can access the list of the providers in this dictionary using the command cs.providers.keys.

# change basemap - can be one of these

# first see available provider names

cx.providers.keys()

dict_keys(['OpenStreetMap', 'OpenSeaMap', 'OpenPtMap', 'OpenTopoMap', 'OpenRailwayMap', 'OpenFireMap', 'SafeCast', 'Thunderforest', 'OpenMapSurfer', 'Hydda', 'MapBox', 'Stamen', 'Esri', 'OpenWeatherMap', 'HERE', 'FreeMapSK', 'MtbMap', 'CartoDB', 'HikeBike', 'BasemapAT', 'nlmaps', 'NASAGIBS', 'NLS', 'JusticeMap', 'Wikimedia', 'GeoportailFrance', 'OneMapSG'])

Once you have the list of providers, you can find the names of their specific tilesets.

Below, we retrieve the list of the tilesets available from the provider CartoDB.

# Get names of tile sets for a specific provider

cx.providers.CartoDB.keys()

dict_keys(['Positron', 'PositronNoLabels', 'PositronOnlyLabels', 'DarkMatter', 'DarkMatterNoLabels', 'DarkMatterOnlyLabels', 'Voyager', 'VoyagerNoLabels', 'VoyagerOnlyLabels', 'VoyagerLabelsUnder'])

Now we can specify a different tileset using the source argument to the add_basemap method.

# Change the basemap provider and tileset

ax = bart_gdf.to_crs('EPSG:3857').plot(figsize=(9, 9))

cx.add_basemap(ax,  source=cx.providers.CartoDB.Positron)

LEARNING MORE

Above, we prove a very short introduction to the excellent contextily library. You can find more detailed information on the contextily homepage, available at: https://github.com/geopandas/contextily. We especially encourage you to check out the notebook examples provided in the contextily github repo. Also, a big shoutout to D-Lab friend Chris Holdgraf in the Division of Computing, Data Science and Society who has been a contributer to this Python library.