Creating bike routes with python

This weekend my goal was to ride 50 miles to and from my house. In my last post I showed four ways to find where I could get to for a certain distance, but I really hate “there and back” rides, so I wanted to find loop-based routes that would have my target distance and not have any doubling back. I used basically the same tools and I’m decently happy with the results.

tl;dr here’s the route

How’d I do it?

Here’s the basic gist of what I did:

Get the lat/long of my house to have a starting point
Do some geo-based math to find points on a regular polygon that includes my house. I set the polygon perimeter length to my target length
Ask openrouteservice to find directions to all of the points in order (beginning and ending at my house)
Adjust the total polygon perimeter until the actual path is my target distance
Construct a google maps URL so that I can navigate while riding

Step one: lat long of my house

This one is easy. Just open google maps and right click anywhere to get the latitude and longitude

After doing that here’s what’s in your clipboard: 44.9240104020827, -93.09348080561426

So just paste that somewhere in your python code as a list and you’re good to go.

Step 2: Geo-based math for other points

I’m pretty sure I could have used some package to do this, but I figured it wouldn’t be that hard to just code up a sloppy version, especially since I knew I’d have to make adjustments to the polygon positions so that the actual travel distance would be what I wanted (see step 4).

If you have the lat/long of a point and want to find a point some known distance away in some particular direction (say $\theta$ degrees counterclockwise from east), you need to realize that while moving north/south follows a great circle (with a radius equal to the earth’s), moving east/west is not moving along a great circle (which is why you have to turn slightly north when driving straight east). Instead you’re moving on a circle whose radius is the earth’s radius multiplied by the cosine (ugh) of the latitude. If you were thinking it should have been sine (like my first thought) it’s likely because you’re used to using polar angles in mathematics (ie, the zero of latitude is the equator, not the north pole). Also, as you’ll note looking at the lat/long I pasted in above, since I’m roughly at 45 degrees latitude, it doesn’t really matter.

So here’s how I coded it:

earthRad=6378.1
def newPoint(start, angle, distance):
  [long,lat]=start
  horizScale=earthRad*math.cos(lat*math.pi/180)
  newLong= long+distance*math.cos(angle)/horizScale*180/math.pi
  newLat=lat+distance*math.sin(angle)/earthRad*180/math.pi
  return [newLong,newLat]

Ok, so now I just need to make a bunch of points on a regular polygon:

def polygon(sides,start,totalDistance,initialAngle=0):
  theList=[start]
  curLoc=start
  for x in range(sides):
    curLoc=newPoint(curLoc,x/sides*2*math.pi+initialAngle,totalDistance/sides)
    theList.append(curLoc)
  return theList

You might wonder about that “initialAngle” part. I found that being able to tilt the initial polygon helped because I have a big river just east of me that the mapping software doesn’t like (only if you try to do bicycle navigating right into the river).

Step 3: Use openrouteservice to map it out

Here’s the code I used:

import openrouteservice as ors
from openrouteservice import convert
import folium
import math
token="GET YOURS FROM OPENROUTESERVICE"
client=ors.Client(key=token)
home=[PASTE THIS IN FROM GOOGLE MAPS]
directions=client.directions(profile='cycling-regular', coordinates=polygon(7,home[::-1],60,3*math.pi/4))
directions['routes'][0]['summary']['distance']

You can get your own token by signing up (for free!) at https://openrouteservice.org/ It has some limitations regarding how often you can use it, but they’re very generous, at least for this kind of thing.

The last line is what you have to do to get the total distance for the journey. You can see that the result of the “client.directions” command is a highly structured object (it took me a while to figure out how to extract that total distance). You get the map I generated above with this:

map=folium.Map(width=500, height=500,location=home)
folium.GeoJson(convert.decode_polyline(directions['routes'][0]['geometry'])).add_to(map)
map

Now you see why I had to import folium up above.

Step 4: repeat and adjust until you get the distance you want

You’ll note in the code above that I submitted points on a polygon that would have a total perimeter of 60 kilometers, but it returned a path that’s just over 80 km, or 50 miles (my original goal). The reason is mostly due to the taxi-cab geometry I talk about in my 2nd approximation in this post. In other words, there’s no way there are dead straight roads among all the polygon vertices, so you’re going to travel further.

Step 5: create a google maps URL

While I like the programmatic capabilities of openrouteservice, I *love* google maps as a navigation aid while riding. So I wanted to figure out how to get my newly found route into google maps. Solution: use their url api! On that page you can learn how to create URLs with starting points, ending points, way points (points along the journey), and type of directions (I wanted bicycle directions). Here’s how I built that:

pgon=polygon(7,home[::-1],60,1*math.pi/4)
baseUrl="https://www.google.com/maps/dir/?api=1&"
start="origin="+str(home[0])+","+str(home[1])
travelmode="travelmode=bicycling"
waypoints="waypoints="+"|".join([str(x[1])+","+str(x[0]) for x in pgon])
destination="destination="+str(home[0])+","+str(home[1])
baseUrl+start+"&"+travelmode+"&"+waypoints+"&"+destination

That gives you a URL that produces this:

And I’m finally ready to ride!

How it went (the actual ride)

It was a beautiful late October day and I really did have a blast. The mapping worked really well but there are some issues I had to deal with.

If you look at the image above, you’ll probably note the little spur on the eastern edge of the route. Essentially google’s mapping algorithm (really both algorithms if you go back and look at the first image on this post) didn’t know of a way to get to and from that particular waypoint without a little bit of doubling back. I didn’t want to do that, so I just figured I could delete that waypoint (google makes that pretty easy to do). However, if you delete it before you begin, google will try to find a faster way from the prior waypoint to the next one, so I just waited until I got close before deleting it. That worked well.

The other waypoints weren’t so far off a normal path, but I did do a lot of weird extra blocks to get to the exact point I asked for. I probably could have saved a few miles not doing that, but that wasn’t really the point.

Overall I think it worked out pretty well. I’m excited to do it again. I’ll probably just tilt the 7-sided polygon over just a little further to get a completely different ride.

Your thoughts?

Got any feedback for me? Here are some starters for you:

I love this! Can you do it for . . .?
I hate this! It’s obvious you could do this in a much easier fashion by . . .
What do you mean you have to turn a little north to drive straight east?
I think we should redo latitude so that the zero is at the north pole
Why python, I thought you loved Mathematica?
Why are the google and openrouteservice directions slightly different?
Why only 50 miles?
How long did it take you? (8 hours but that was with a really long lunch and a few other stops to read and watch rugby)
How did you find rugby to watch?
What’s wrong with using cosine (you said “ugh” after doing so)?