Fork me on GitHub

Using Folium to show geographic data

This post demonstrates how to

  • Display a map from OpenStreetMaps using Folium
  • Add custom shape files to define regions of interest
  • Color the regions of interest based on data in a pandas dataframe
  • (New) Brief introduction to GeoPandas.

The code of this post is available at https://github.com/rsandstroem/IPythonNotebooks/blob/master/GeoMapsFoliumDemo/GeoMapsFoliumDemo.ipynb .

August 1, 2017: This post is almost a year old, but I decided to make some technical improvements to it to improve the viewing experience online. (It is the Swiss National Day after all!) During this renovation, I stumbled upon GeoPandas, which is a great package worth mentioning in this context.

Getting map data and file formats

Using the method described here, you do not need any local map data. The maps are automatically pulled from OpenStreetMap.

However, you might want use concepts such as regions, and you need to define them through files containing those geolocations. There are multiple sources of data which is freely available, in this example we are using shape files from this location: http://data.geo.admin.ch.s3.amazonaws.com/ch.swisstopo.swissboundaries3d-kanton-flaeche.fill/data.zip

We have provided the shape files we will use for this demonstration, normally you would need to retrieve and unpack the map data.

There are many file formats commonly used for maps. In this case what we have are a set of shapefiles. Since we will be working with geojson files in this demonstration you need to convert them to geojson.

A sample geojson file has been supplied with this tutorial, but for reference you can create a geojson file from the shape files from a console like this:

ogr2ogr -f GeoJSON -t_srs EPSG:4326 -simplify 1000 switzerland.geojson swissBOUNDARIES3D_1_2_TLM_KANTONSGEBIET.shp

You now have a map in your local directory, formatted as a geojson file. I use the simplify option to reduce the file size, at the expense of accuracy.

Set up

For this demonstration we need "folium". Install it using

pip install folium

or if you are using the Anaconda you can do

sudo conda install --channel https://conda.binstar.org/IOOS folium

where you can replace the channel with your choice from the options listed by

binstar search -t conda folium
In [1]:
import folium
print folium.__version__
0.3.0

Display the map

First we show the map without any use of the geojson we created from the shapefiles.

In [2]:
kanton_map = folium.Map(location=[46.8, 8.33],
                   tiles='Mapbox Bright', zoom_start=7)
kanton_map
Out[2]:

You should now see an map of Switzerland with surrounding countries and cities. Next, display an OpenStreetMap and overlay the cantons of Switzerland on top from the geojson file created earlier.

In [3]:
kanton_map.choropleth(geo_path='switzerland.geojson')
kanton_map
Out[3]:

Import demographics data from a csv into a pandas data frame

The purpose of using a color filled map like the one above is usually to have the color represent something in the data, e.g., sales revenue per district. Since this data is separate from the geometrical description we need to retrieve it from an external source and join it with the geometrical shapes from the geojson file.

Read in some data for the canton's of Switzerland from a csv file. We will display this data on the map we just created.

In [4]:
import pandas as pd
kanton_overview = r'Switzerland_overview.csv'
kanton_data = pd.read_csv(kanton_overview)
kanton_data
Out[4]:
Abbr Canton Capital Population Area Density Municipalities Official languages CantonNumber
0 ZH Zurich Zurich 1463459 1729 701 171 German 1
1 BE Bern Bern 1001281 5959 158 383 German, French 2
2 LU Luzern Luzern 390349 1493 233 87 German 3
3 UR Uri Altdorf 35865 1077 33 20 German 4
4 SZ Schwyz Schwyz 151396 908 143 30 German 5
5 OW Obwalden Sarnen 36507 491 66 7 German 6
6 NW Nidwalden Stans 41888 276 138 11 German 7
7 GL Glarus Glarus 39593 685 51 3 German 8
8 ZG Zug Zug 118118 239 416 11 German 9
9 FR Fribourg Fribourg 297622 1671 141 167 French, German 10
10 SO Solothurn Solothurn 261437 790 308 122 German 11
11 BS Basel-Stadt Basel 195962 37 5072 3 German 12
12 BL Basel-Landschaft Liestal 281112 518 502 86 German 13
13 SH Schaffhausen Schaffhausen 78783 298 246 27 German 14
14 AR Appenzell Ausserrhoden Herisau 53691 243 220 20 German 15
15 AI Appenzell Innerrhoden Appenzell 15778 173 87 6 German 16
16 SG St. Gallen St. Gallen 491699 2026 222 85 German 17
17 GR Graubünden Chur 194959 7105 26 180 German, Romansh, Italian 18
18 AG Aargau Aarau 636362 1404 388 220 German 19
19 TG Thurgau Frauenfeld 260278 991 229 80 German 20
20 TI Ticino Bellinzona 346539 2812 110 157 Italian 21
21 VD Vaud Lausanne 749373 3212 188 339 French 22
22 VS Valais Sion 327011 5224 53 143 French, German 23
23 NE Neuchatel Neuchatel 177327 803 206 53 French 24
24 GE Geneva Geneva 479158 282 1442 45 French 25
25 JU Jura Delemont 72410 839 82 64 French 26

Color code the cantons based on their population density

The KANTONSNUM in the json file is used to ensure correct assignment of the data in the demographics data file. Try using the name of the cantons instead. What happens, and why? (Hint: open the geojson file and compare with the csv file.)

In [5]:
kanton_map2 = folium.Map(location=[46.8, 8.33], 
                    zoom_start=7.5)
kanton_map2.choropleth(geo_path='switzerland.geojson', data=kanton_data,
             columns=['CantonNumber', 'Density'],
             key_on='feature.properties.KANTONSNUM',
             threshold_scale=[100, 200, 300, 500, 1000, 5000],
             fill_color='BuPu')
kanton_map2
Out[5]:

Using GeoPandas

While making some technical improvements to this post, I stumbled upon GeoPandas. While I cannot say I have used it much yet, I think it is a great way to read in geojson files and work with their contents.

In [6]:
import geopandas as gpd
df = gpd.read_file('switzerland.geojson')
print df.columns
df.head()
Index([u'DATUM_AEND', u'DATUM_ERST', u'EINWOHNERZ',  u'ERSTELL_J',
        u'ERSTELL_M', u'GRUND_AEND',   u'HERKUNFT', u'HERKUNFT_J',
       u'HERKUNFT_M',        u'ICC', u'KANTONSFLA', u'KANTONSNUM',
          u'KT_TEIL',       u'NAME',  u'OBJEKTART', u'REVISION_J',
       u'REVISION_M', u'REVISION_Q', u'SEE_FLAECH',       u'UUID',
         u'geometry'],
      dtype='object')
Out[6]:
DATUM_AEND DATUM_ERST EINWOHNERZ ERSTELL_J ERSTELL_M GRUND_AEND HERKUNFT HERKUNFT_J HERKUNFT_M ICC ... KANTONSNUM KT_TEIL NAME OBJEKTART REVISION_J REVISION_M REVISION_Q SEE_FLAECH UUID geometry
0 2014-10-31 2012-10-26 194935 2012 10 verbessert swisstopo 2014 8 CH ... 18 0 Graubünden Kanton 2015 1 1412_Aufbau 0.0 {0B2364ED-49E0-4D53-A33C-C684DD530B57} POLYGON Z ((8.878214669955852 46.8142107947290...
1 2014-10-31 2012-10-26 1000918 2012 10 verbessert swisstopo 2014 8 CH ... 2 1 Bern Kanton 2015 1 1412_Aufbau 11897.0 {DDD56CEF-0E61-4EED-85ED-F67A459C93ED} POLYGON Z ((7.237992928187036 46.5550699033581...
2 2014-10-31 2012-10-26 326553 2012 10 verbessert swisstopo 2014 8 CH ... 23 0 Valais Kanton 2015 1 1412_Aufbau 1060.0 {54B25E50-30A7-4995-ADE3-5FFF6E13A995} POLYGON Z ((8.411360233490436 46.6542191237368...
3 2014-10-31 2012-10-26 748815 2012 10 verbessert swisstopo 2014 8 CH ... 22 1 Vaud Kanton 2015 1 1412_Aufbau 39089.0 {921DFEF2-6D91-4CB8-9CFC-2A831C412020} POLYGON Z ((6.780677661693769 46.8542740178914...
4 2014-10-31 2012-10-26 346472 2012 10 verbessert swisstopo 2014 8 CH ... 21 0 Ticino Kanton 2015 1 1412_Aufbau 7152.0 {95F10F52-8B2F-4D6A-AF7E-D4F915E42F89} POLYGON Z ((8.478723733515089 46.5288928701094...

5 rows × 21 columns

Very easy to make a pandas dataframe from this json file! Every row represents one area, and the last column contains the shape of the area as a three dimensional polygon, i.e., the longitude, latitude and altitude of the regional borders. The dataframe also contains data columns, such as number of inhabitants (EINWOHNERZ) and surface area (KANTONSFLA).

A nice feature of using GeoPandas in a Jupyter Notebook is the ease at which we can draw the content of the dataframe:

In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12,8))
ax.set_aspect('equal')
df = df.to_crs({'init': 'epsg:4326'}) # change projection
df['density'] = df.EINWOHNERZ / df.KANTONSFLA # calculate population density
df.plot(ax=ax, column='density', cmap='OrRd')
Out[7]:

Since the geometry of the cantons are represented as Shapely polygons, Jupyter can render them directly. For example, the canton of Valais has index 2 in the dataframe df:

In [8]:
print type(df.geometry[2])
df.geometry[2]

Out[8]:

Example: Find mountains along the borders.

Let's combined what we learned and add map pins for mountains on the border where the altitude exceeds 3500 m. First we extract the three dimensional points from the dataframe where the altitude exceeds 3500 m.

In [9]:
import numpy as np

peaks = list()
for kanton in df.index:
    arr = np.array(df.geometry[kanton].exterior.coords)
    peaks.append(arr[np.where(arr[:,2]>3500)].tolist())
peaks = [item for sublist in peaks for item in sublist] # flatten to a single list of 3d points 
print peaks[:3]
[[8.91603243665709, 46.812444938780914, 3604.6812499999796], [9.960072829925547, 46.37943103812307, 3901.2637499999837], [9.935780191298486, 46.37607724456639, 3889.3874999999825]]

Then we create a feature group consisting of map markers for the peaks we found in the previous step. The altitude is added as a pop-up that will be displayed if the user clicks on a map marker. Finally we add this new feature group to the populations density map which we created earlier.

In [10]:
feature_group = folium.FeatureGroup("Locations")
for peak in peaks:
    feature_group.add_child(folium.Marker(location=[peak[1], peak[0]], popup=str(peak[2]))) 
kanton_map2.add_child(feature_group)
kanton_map2
Out[10]:

The regions with high border mountains are located along the southern border of Switzerland, the Alpes. According to the population density we overlayed earlier, the cantons with high mountains are also the least populated. Hardly surprising, it is hard to commute to work if you live on top of a glacier!

Share on: LinkedInTwitterFacebookGoogle+Email

Comments !

blogroll

social