Choropleth Maps in Python (2021)

First maps with geoviews¶

We import the tools we need, and we initialize geoviews:

geoviews.feature contains the geometric description of various features, such as the world land, ocean, and borders. Let's plot this:

With cartopy, we can easily change the projection:

Finding more features¶

For now, we have extracted the features we wanted to plot (gf.land, gf.ocean, gf.borders) from geoviews.features. Let's see what else we can find there:

In fact, these are only pre-defined features for common use cases, and are quite limited. For example, the states features are restricted to the US.

Cartopy, which is in fact used under the hood by geoviews to get the land, ocean, and borders features, has a feature interface that allows to get access to more features through a variety of data sources, like (Natural Earth)(http://www.naturalearthdata.com).

But if you're looking for something specific, it might be unavailable. For example, I haven't been able to find a geographic description of the french departments, to create the map at the top of this article.

In the next section, we will see how to work around this problem.

Loading a GeoJSON file with pandas¶

GeoJSON files encode geographic data structures and can easily be read and manipulated with geopandas.

A google search of "france geojson" led me to the france-geojson repository by Grégoire David (thanks!). From this repository, I downloaded departements-version-simplifiee.geojson as france-geojson/departements-version-simplifiee.geojson.

After doing that, we can import geopandas and read this file:

We see that the department geometries are stored as polygons, with the longitude and latitude of each point on the polygon. It's easy to draw a department with geopandas: we just select a row from the dataframe (here the first one for the department Ain), and we evaluate the geometry:

We're now going to show all french departments, with a color indicating a measured value for each department. For the moment, we don't have a measured value to display, so we generate one randomly:

Let's have a closer look at this line of code.

On the right side, we use numpy to generate an array of random integers between 1 and 10, with a length equal to the number of rows in the dataframe, which is sf.shape[0]. Then, we assign this numpy array to a new column in the dataframe, called value.

Plotting geographical data with geoviews¶

In the holoviews style, we just need two lines to display the departments. The first line declares the set of polygons from our geopandas dataframe. We explicitely choose the department name and its value as value dimensions:

Then, we plot the departments, with a few options. The most important one is color=dim('value'). This just means that we are going to map value to a linear color scale for each department. And you can find a full explanation of what dim objects are in the holoviews documentation.

Note that the column used for the color must be in vdims. It's easy to forget!

We used bokeh as a backend, so we benefit from the usual tools from this library (at the top of the plot, as specified), and from its interactive features. You can try to wheel zoom, pan, and hover on the departments to get detailed information.

Creating a choropleth map with your own data¶

In the simple example above, we randomly defined a value for each department. But in real life, what you want to do is to display data that is not present in the geopandas dataframe. So we need a way to associate data to the geopandas dataframe before plotting.

That's what we're going to do now, taking as an example the estimated population in each department, in 2019. I found the data here in xls format, and I converted it to tsv file (tab-separated values) to create the file population_france_departements_2019.tsv.

To read this file, we use plain pandas. Indeed, we don't need geoviews here because this dataset does not contain any geographical information. We use the well-known read_csv method, specifying that tabs are used to separate values, and that a comma is used to indicate thousands. You can just have a look inside the tsf text file if you want.

After the first columns that indicate the deparment code and name, we have three sets of columns for the number of people in the department. The prefixes mean the following:

• t_ : total number
• h_ : number of men
• f_ : number of women

And in each set, there is one column for each age range: from 0 to 19 years, ... , over 75 years.

Now, we need to associate the data in this pandas dataframe (df) to our geopandas dataframe (sf). That's very easy with the dataframe merge method.

Below, we merge data on each line, making sure that the value in column code of sf is equal to the value in column depcode of df. The suffixes option is for columns that have the same name in df and sf, like nom (the department name). It determines the names that the columns will have after merging. Here, I want the column that came from sf to be kept with the same name, and the one that came from df to get a _y suffix.

And we're now ready to plot the total population in each department! We also add the total number of men and women in the tooltips, by adding these columns to vdims. Also, please note how I rescaled the t_total color scale to display it in millions of inhabitants.