How to find places along a GPX route in Python using Overpass API and Google Maps Places API

Marrying free and paid APIs

Author

Andrew Wang

Do you have some long GPX routes, for long hikes, pub crawls or pilgrimages? How can you automatically find what amenities are near the route, and reduce the amount you spend using paid APIs?

In this tutorial, we find places in OpenStreetMap along our route, using the free Overpass API. Then, we use the Google Maps Places API for each place of interest to retrieve up-to-date place ratings, reviews, photos and opening hours.

Doing it this way has a couple of advantages:

What we need to do is:

  1. Import a GPX route into Python
  2. Find nearby places (for free, using Overpass API)
  3. Find details of places (using Google Maps Places API)
  4. Display the places on a route map!

For the full notebook, see the repo.

Let’s get started!

Create a virtual environment: python -m venv venv, source venv/Scripts/activate

You will need to pip install the following packages (pip install -r requirements.txt):

import numpy as np
import pandas as pd

1. Import GPX route

We will be using a sample GPX route describing a pilgrimage in the Peak District, United Kingdom, which you can download from the British Pilgrimage Trust website. The route is 63km long and wanders through small villages and beautiful countryside.

After downloading the route, read the GPX using the gpxpy package and convert to a pandas DataFrame (this handles both routes and tracks):

from gpxpy import parse

with open("Peak Pilgrimage.gpx") as f:
    p = parse(f)
    points = [(point.latitude, point.longitude) for route in p.routes for point in route.points] + \
             [(point.latitude, point.longitude) for track in p.tracks for segment in track.segments for point in segment.points]
    input_df = pd.DataFrame(points, columns=["latitude", "longitude"])

The route might look like this:

input_df.head()
latitude longitude
0 53.053368 -1.803243
1 53.053603 -1.803120
2 53.053691 -1.802744
3 53.053987 -1.802573
4 53.054065 -1.802433

Plot the route on a map:

Python function to plot a GPX route and POIs to a HTML map using Folium 
import folium

def folium_route(route_lats: pd.DataFrame, 
                 route_lons: pd.DataFrame, 
                 pois=[], popup_func=lambda x,y: None) -> folium.folium.Map:
    
    # Initialise map at centre of route
    m = folium.Map([route_lats.mean(), route_lons.mean()], zoom_start=10)
    
    # Add route to map
    folium.PolyLine(list(zip(route_lats, route_lons))).add_to(m)
    
    # If displaying POIs, plot as circles with popups
    for i in range(len(pois)):
        folium.CircleMarker(location=pois.loc[i, ["lat", "lon"]].to_list(), fill=True) \
              .add_child(folium.Popup(popup_func(pois, i), max_width=500))             \
              .add_to(m)
        
    return m
folium_route(input_df["latitude"], input_df["longitude"])
Make this Notebook Trusted to load map: File -> Trust Notebook

2. Find nearby places along the route using Overpass API

We use the free Overpass API to find places from OpenStreetMap along our GPX route. The overpy library lets us do this in Python:

import overpy 
api = overpy.Overpass()

Flatten the route coordinates into a string:

input_latlons = ",".join(input_df.to_numpy().flatten().astype("str"))

Construct the Overpass query in Overpass Query Language. See here for a great intro.

  • nwr tells Overpass to find nodes, ways and relations (explanation).
  • amenity is a key-value tag describing facilities such as pubs, see here for a complete list.
  • historic describes features of historic interest, see here for a complete list.
  • around:1000 finds features within 1000m of the input route.
overpass_query = f"""
[out:json][timeout: 500];
(nwr["amenity" = "place_of_worship"](around:1000,{input_latlons});
nwr["amenity" = "pub"](around:1000,{input_latlons});
nwr["historic" = "memorial"](around:1000,{input_latlons});
nwr["historic" = "wayside_cross"](around:1000,{input_latlons});
nwr["historic" = "wayside_chapel"](around:1000,{input_latlons});
nwr["historic" = "wayside_shrine"](around:1000,{input_latlons});
nwr["historic" = "archaeological_site"](around:1000,{input_latlons});
nwr["historic" = "stone"](around:1000,{input_latlons});
nwr["historic" = "milestone"](around:1000,{input_latlons});
nwr["historic" = "boundary_stone"](around:1000,{input_latlons});
nwr["historic" = "tomb"](around:1000,{input_latlons});
nwr["historic" = "archaeological_site"](around:1000,{input_latlons});
nwr["water_source" = "well"](around:1000,{input_latlons});
  ( ._; >; );
);
out center;
"""
result = api.query(overpass_query)

The Overpass result is returned as a list of nested objects, which we need to map to a json dictionary. The following pandas steps then loads the results into a DataFrame and removes all results with an empty name field:

result_df = pd.json_normalize(map(vars, result.nodes)) \
              .drop(columns=["_result"]) \
              .rename(columns={"tags.name": "name"}) \
              .dropna(subset=['name']) \
              .reset_index(drop=True)

We see that Overpass gives us a rich set of results directly from OpenStreetMap:

result_df.columns
Index(['id', 'lat', 'lon', 'tags.addr:city', 'tags.addr:postcode',
       'tags.addr:street', 'tags.addr:village', 'tags.amenity', 'tags.fhrs:id',
       'name', 'tags.food', 'tags.outdoor_area', 'tags.parking',
       'tags.real_ale', 'tags.source:postcode', 'tags.wikidata',
       'tags.accommodation', 'tags.floor:material', 'tags.internet_access',
       'tags.internet_access:fee', 'tags.toilets', 'tags.toilets:wheelchair',
       'tags.website', 'tags.wheelchair', 'tags.opening_hours:signed',
       'tags.operator', 'tags.note', 'tags.archaeological_site',
       'tags.historic', 'tags.site_type', 'tags.source', 'tags.wikipedia',
       'tags.denomination', 'tags.ele', 'tags.religion', 'tags.memorial',
       'tags.tourism', 'tags.megalith_type', 'tags.addr:housename',
       'tags.sport', 'tags.inscription', 'tags.heritage'],
      dtype='object')
result_df.head()
id lat lon tags.addr:city tags.addr:postcode tags.addr:street tags.addr:village tags.amenity tags.fhrs:id name ... tags.denomination tags.ele tags.religion tags.memorial tags.tourism tags.megalith_type tags.addr:housename tags.sport tags.inscription tags.heritage
0 251982730 53.2459014 -1.6133234 Bakewell DE45 1SR Nether End Baslow pub 258157 The Devonshire Arms ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1 251985589 53.2465269 -1.6124122 Bakewell DE45 1SR Nether End Baslow pub 258453 The Wheatsheaf ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 251986045 53.2661742 -1.6310557 NaN S32 3XA Calver Bridge NaN pub 258965 The Bridge Inn ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
3 251986258 53.2683640 -1.6411023 NaN NaN NaN NaN pub NaN The Derwentwater Arms ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
4 305916852 53.2814117 -1.6312423 NaN NaN NaN NaN pub NaN The Chequers ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

5 rows × 42 columns

Click on the points of interest (POI) to see their name:

folium_route(input_df["latitude"],
             input_df["longitude"], 
             pois=result_df, 
             popup_func=lambda p,i: f"POI {i}: {p.loc[i,'name']}")
Make this Notebook Trusted to load map: File -> Trust Notebook

The eagle-eyed would notice that the POIs aren’t listed in the order of the route! This can present some problems. Let us use reorder them using some distance-based matching:

  1. Calculate pairwise distances between POIs and route points
  2. Find closest route point for each POI
  3. Reorder POIs by closest route points
from sklearn.metrics.pairwise import haversine_distances as hd

pairwise_distances = hd(np.deg2rad(input_df[["latitude", "longitude"]].to_numpy()), 
                        np.deg2rad(result_df[["lat", "lon"]].to_numpy())) \

result_df["closest"] = pairwise_distances.argmin(axis=0)

pois = result_df.sort_values("closest")
pois = pois[["id", "lat", "lon", "name"]].reset_index(drop=True)

3. Get full details of POIs using Google Maps Places API

For each of our POIs on the route, we want to quickly grab their website, photos, and ratings. Inevitably, OpenStreetMap has limited and incomplete data on all of our detected places. Here we use the Google Maps Places API to get those details.

To get started with the API, grab an API key from the Google Maps Platform (paid service). Note that all projects get a free monthly allowance which will be more than enough for us: see billing for details.

The googlemaps Python SDK allows us to connect to the Places API in Python (copy-paste your personal API key in):

import googlemaps

gmaps = googlemaps.Client(key="")

First we perform a “Nearby Places” search to get the corresponding Google Maps place for each OSM place, by providing Google Maps with both the location and the OSM name of the place.

poi_lat_lon_names = list(pois[["lat", "lon", "name"]].itertuples(index=False, name=None))
from tqdm import tqdm

results = [gmaps.places_nearby(keyword=name, location=[lat, lon], radius=5) \
           for (lat, lon, name) in tqdm(poi_lat_lon_names)]
100%|███████████████████████████████████████████| 26/26 [00:14<00:00,  1.75it/s]

Coax the results into a clean json format for loading into a DataFrame:

results = [result["results"][0] if len(result["results"]) > 0 else {"business_status": np.nan} for result in results]
gmaps_df = pd.json_normalize(results)

Now, extract useful fields for our POIs:

  • Name on Google Maps
  • Average user rating and number of rating
  • Type of place according to Google Maps
  • Price level
  • Google URL
  • User uploaded photo
pois[["google_name", "rating", "types", "n_ratings", "price_level"]] = \
    gmaps_df[["name", "rating", "types", "user_ratings_total", "price_level"]]

pois["url"] = gmaps_df["place_id"].map(lambda x: f"https://www.google.com/maps/place/?q=place_id:{x}")

Perform a “Place Photos” search to obtain user photos from Google Maps. The photos are returned in byte chunks.

from io import BytesIO
import base64

def get_photo_bytes(photo: list) -> str:
    if isinstance(x, list) and len(x) > 0:
        
        # Get top returned photo ID
        photo_id = photo[0]["photo_reference"]
        
        # Get raw bytes of photo from Google Maps Place Photos search
        b = b"".join(gmaps.places_photo(photo_id, max_width=100))
        
        # Encode into base64 for HTML display
        return base64.b64encode(BytesIO(b).getvalue()).decode()        
    else:
        return ""
pois["photo_bytes"] = gmaps_df["photos"].map(get_photo_bytes)

Google Maps has given us a lot more detail than the original OSM places!

pois.head()
id lat lon name google_name rating types n_ratings price_level url photos photo_id photo_bytes
0 3362045361 53.141174 -1.809284 Royal British Legion Royal British Legion, Cavendish House 4.8 [bar, point_of_interest, establishment] 8.0 NaN https://www.google.com/maps/place/?q=place_id:... AW30NDzDoXCGbKue-GlnzehJqq02gCrxEeV5oBn0rrZYBq... AW30NDzDoXCGbKue-GlnzehJqq02gCrxEeV5oBn0rrZYBq... /9j/4AAQSkZJRgABAQAAAQABAAD/4QAqRXhpZgAASUkqAA...
1 410871346 53.140523 -1.809675 Devonshire Arms The Devonshire Arms 4.1 [bar, restaurant, food, point_of_interest, est... 627.0 2.0 https://www.google.com/maps/place/?q=place_id:... AW30NDxJj-kzaTz5yIR9sfw9mFd3ApawYQynbViq0VNPLE... AW30NDxJj-kzaTz5yIR9sfw9mFd3ApawYQynbViq0VNPLE... /9j/4AAQSkZJRgABAQAAAQABAAD/4QAqRXhpZgAASUkqAA...
2 3362042956 53.140606 -1.808380 War Memorial NaN NaN NaN NaN NaN https://www.google.com/maps/place/?q=place_id:nan NaN NaN
3 8289562776 53.193764 -1.826558 Cronkston Low Bowl Barrow NaN NaN NaN NaN NaN https://www.google.com/maps/place/?q=place_id:nan NaN NaN
4 414090489 53.195983 -1.776546 The Bulls Head The Bulls Head 4.6 [bar, restaurant, point_of_interest, food, est... 1082.0 2.0 https://www.google.com/maps/place/?q=place_id:... AW30NDxZpxIryMv0BWcYDKa-SIe1_LYizfasNIMfgwh3-2... AW30NDxZpxIryMv0BWcYDKa-SIe1_LYizfasNIMfgwh3-2... /9j/4AAQSkZJRgABAQAAAQABAAD/4QAqRXhpZgAASUkqAA...

4. Display final map

Time to display the final results in a map. Click on the POIs to get a full popup with Google Maps details and photos! Note some POIs weren’t detected in Google. These are, for example, war memorials that a user has added to OSM but haven’t yet been added to Google!

Python function to return pretty HTML popups for Folium map 
def popup_func(pois: pd.DataFrame, i: int) -> str:
    
    if not isinstance(pois.loc[i, 'google_name'], str):

        # POI with no Google results
        return f"POI {str(i)}: {pois.loc[i, 'name']}"
    else:
        
        # Pretty HTML for POI details
        return f"""
POI {str(i)}: <a href=\"{pois.loc[i, 'url']}\">{pois.loc[i, 'google_name']}</a><br>
Rating: {pois.loc[i, 'rating']}<br>
Number of ratings: {int(pois.loc[i, 'n_ratings'])}<br>
Tags: {', '.join(eval(pois.loc[i, 'types'])).replace('_', ' ')}<br>
<img src=\"data:image/jpeg;base64,{pois.loc[i, 'photo_bytes']}\">
    """
folium_route(input_df["latitude"], 
             input_df["longitude"], 
             pois=pois,
             popup_func=popup_func)
Make this Notebook Trusted to load map: File -> Trust Notebook