Gemello Engineering AB - Compute This!

This time of the year the new year resolutions often still apply. For some people these resolutions often involve striving for better choices. This can be working out at a gym (if the Corona restrictions allow this) or to make more sensible choices when it comes to what dinner to eat. This blog post will be on the opposite, i.e. how to get hold of some good old fashioned junk food.

Where to find the golden arches?

Growing up in “the countryside” in the southern part of Sweden (Skåne), there was no such thing as comercial junk food in the close by areas (not if you were riding a bike). Since I have been a sucker for junk food my entire life, this was of course quite annoying. When the craving for a Big Mac was present, I felt like we lived in the most junk food hostile area in all of Skåne. But was this true? This leads to the subject for this post:

Does there exist a place in Skåne where the distance to junk food, in this case McDonald's, is greater than the place I grew up at?

To try to find an answer to this question, let's fire up a Jupyter notebook and get ready to do some coding.

How to work with maps?

Python has a few ways to work with maps. In this post we will use a number of different packages designed to do this. Let's start with doing some imports:

import os
import numpy as np
from geopy.distance import geodesic 
import pandas as pd

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import maskoceans

Now let’s plot a map of Skåne and mark out some cities:

# Load data for city
City = pd.read_csv('data/City.csv',sep=';')

fig = plt.figure(figsize=(8, 8))
m = Basemap(
        projection='merc',
        resolution='h', 
        area_thresh=250,
        llcrnrlon=12.0, 
        llcrnrlat=55.3,
        urcrnrlon=15.0, 
        urcrnrlat=56.5)
m.drawcoastlines()
m.drawmapboundary(fill_color='lightgray')
m.fillcontinents(color='white',lake_color='aqua')

for i in range(len(City)):
    x, y = m(City.iloc[i].lng,City.iloc[i].lat)
    plt.plot(x,y,'bo',alpha=0.5)
    plt.text(x+3000, y,City.iloc[i].Name,fontsize=10)

This will result in the following map:

That looks really nice. Here I have also added the place where I grow up (Klamby) which I can assure you is not a city…

Let's get some data…

OK, the next step will now be to get hold of data for locations of McDonald's in Skåne. So how to do this? One approach could be to go online and do a manual search using some online map services. However, we will try some other approach here and use the API towards Open Street Map. In particular, we will use the Overpass API in combination with python.

Before we start, we have to take a look at how Open Street Map is structured. We have three basic components in the OSM data model, which are nodes, ways and relations which all come with an id. Many of the elements come with tags which describe specific features represented as key-value pairs. In simple terms, nodes are points on the maps (in latitude and longitude). A way on the other hand is a ordered list of nodes, which could correspond to a street or the outline of a house. The final data element is a relation which is also an ordered list containing either nodes, ways or even other relations. It is used to model logical or geographic relationships between objects.

OK, that was some introduction to how Open Street Map is organised, but how do we search for the locations of McDonald's using this interface? Well, looking at the wiki the following tagging is recommended amenity=fast_food so why not start with that? In addition, we will utilize the requests command to perform the quiery. The python code looks something like this:

import requests
import json
overpass_url = "http://overpass-api.de/api/interpreter"
overpass_query = """
[out:json];
(node["amenity"="fast_food"](55.3,12.0,56.5,15.0);
 way["amenity"="fast_food"](55.3,12.0,56.5,15.0);
 rel["amenity"="fast_food"](55.3,12.0,56.5,15.0);
);
out center;
"""
response = requests.get(overpass_url, 
                        params={'data': overpass_query})
data = response.json()

Here we got all the fast food places, and not just the McDonald's places. To get hold of the coordinates for the different locations some additional post processing is needed. This is done below:

# Collect coords into list
kk = 0
coords = []

for element in data['elements']:
    kk +=1
    if element['type'] == 'node':
        lon = element['lon']
        lat = element['lat']
        coords.append((lon, lat))
    elif 'center' in element:
        lon = element['center']['lon']
        lat = element['center']['lat']
        coords.append((lon, lat))

# Convert coordinates into numpy array
X = np.array(coords)

Now to find the places related to McDonald's, we need to do an additional search. Here we simply look thou all of the tags and searching for part of the company name as shown below:

mc_index = []

for element in data['elements']:
    OK = 0
    for ky in element['tags'].keys():
        if 'onald' in element['tags'][ky]:
            OK = 1
    if OK == 1:
        mc_index.append(True)
    else:
        mc_index.append(False)
    
mc_index = np.array(mc_index)    
X_MC = X[mc_index]

Sweet, now we (hopefully) have the locations of the different restaurants having those golden arches on top of the building.

The final step is now to create a grid of locations and then find the closest McDonald's from that location. Let's create a function that will return the distance to the closest restaurant given an arbitrary location.

def GetMacDist(X_MC,Poss):
    close = 100000
    for i in range(len( X_MC)):
        dist = (geodesic(Poss, (X_MC[i][1],X_MC[i][0]) ).km)

        if dist < close:
            close = dist
    
    return close

Now let's do the final steps and just do the hard work of going through all points in the grid and compute the distance to the closest McDonald's (and also do a check if the investigated point is on land or not).

dist_vec = []
size_vec = []
fig = plt.figure(figsize=(8, 8))

# Loop over the grid and find the closest restaurant
for i in range(len(x_vec)):

    if m.is_land(x_vec[i],y_vec[i]):
        temp = m(x_vec[i],y_vec[i],inverse=True)
        dist_vec.append( GetMacDist(X_MC, (temp[1],temp[0]) ) )
        size_vec.append(30)
    
    else:
        dist_vec.append(0)
        size_vec.append(0)

# Do a scatter plot to generate the heat map        
plt.scatter(x_vec,y_vec,s=size_vec,c=dist_vec)

# Draw the cost line
m.drawcoastlines()

And finally lets create a heat map of the distance to the closest McDonald's to see if Klamby is the most junk food hostile area in Skåne or there exist some other point? The heat map shows the distance to the closest McDonald's in km.

Allright, this looks pretty cool. Klamby is quite bad actually. A bit surprisingly, Kivik turns out to be a spot far from the golden arches. But who needs hamburgers, when you can have a Sillaburgare (Heringburger)!!

Time to conclude

So there we have it. The point where I grew up is not the most junk food hostile area in Skåne, but still bad. Good to know! Hope you like this blog post and as always, please feel free to contact us on info@gemello.se.

Till next time, stay safe and have a burger.