Well here we are! Time to solve our dilemma we begun last week! If you haven’t already you should check out Part 1 of this blog which covers the Haversine formula we use for calculating distances. To quickly recap – our problem is that we have a hypothetical app, in our app we have a lot of users, and we want to show the users all of the restaurants with a 30 km radius. Our app needs to work in more than 1 country, so no smart postcode type system is available, and we are expecting to have millions of potential users and restaurants. Basically we don’t want to filter through every restaurant and complete the Haversine formula for each!

So how did we challenge this, our main goal is to move some of the process onto the database so we are only extracting records that either are in our search radius, or are very close to – in our example of the Haversine formula from part 1 we had a user in the middle of Utah, USA. In the case of this user it would be fair to say that we could limit our restaurants to places only within the State of Utah; but in our client’s case we weren’t talking about the USA, and while we still could have looked at regions, let’s travel to the South East corner of Utah, selecting only restaurants in Utah might exclude 270° of our 30 km radius which falls in Colorado, New Mexico and Arizona. So I decided this wasn’t acceptable.

Of course it wouldn’t be too far fetched to say let’s fetch the current state and those around it – but again looking at the USA, there is a point in the far North of West Virginia where the width of the state is less than 15 km, this means that standing on the far east of Ohio, your 30 km circle would include West Virginia AND Pennsylvania, the next state over as it were. Of course we could write exceptions for places like this, but to then cover the world? This is getting out of hand.

Finally I came up with the following solution. The first question is how I chose to store the data in the database – thankfully this is a new project so no data migration was required. I decided to split the longitude and latitude into two different columns. This allowed me to do the following roughly:

Yup that’s the simple solution. Of course this essentially draws a skewed square around the user meaning in the far north west corner you would be more than 30 kms away – but using this method I can substantially reduce the number of potential restaurants to filter through using the Haversine formula knowing that over 80% are going to match! But how do I do all this?

First we need to reverse the Haversine formula. We now know the distance – and the bearing we need to travel in, we just need to find the point!

Keep in mind the earth has a radius of 6,371 km, so to get the radians for 30 kms we need to calculate 30 / 6371. We also need to use our compass degrees and convert them into radians. We will also convert our start point into radians from the equator and the prime meridian. Again we will use our example of Swasey Peak, Utah, USA (39.388281, -113.316093).

// Again we are rounding and only using 5 decimal places
var distance = 30 / 6371;
// distance = 0.00471
var nRad = 0; // North = 0° so very simple
var eRad = 90 * Math.PI / 180;
// eRad = 1.57080
var sRad = 180 * Math.PI / 180;
// sRad = 3.14159 (PI)
var wRad = 270 * Math.PI / 180;
// wRad = 4.71239
var latR = 39.388281 * Math.PI / 180;
// latR = 0.68745
var lonR = -113.316093 * Math.PI / 180;
// lonR = -1.97774

Right – now we are set up with all the radian figures we will need. We have to repeat the process 4 times, one for each direction (North, East, South and West). For this demonstration I am only going to do it for East (As North is too simple since it’s a 0!). To get the latitude of our new point we use the sine of our starting latitude, multiplied by the cosine of our distance, and then add it to the value of the cosine of our latitude multiplied by the sine of our distance multiplied by the cosine of our bearing. So let’s look at that in code for east – bearing in radians of 1.57080. Then take the inverse sine of all of that!

// New latitude function asin(sin(latR) * cos(distance) + cos(latR) * sin(distance) * cos(eRad))

var newLat = Math.asin(Math.sin(latR) * Math.cos(distance) + Math.cos(latR) * Math.sin(distance) * Math.cos(eRad));
// newLat = Math.asin( Math.sin(0.68745) * Math.cos(0.00471) + Math.cos(0.68745) * Math.sin(0.00471) * Math.cos(1.57080) )
// newLat = Math.asin( 0.63457 * 0.99999 + 0.77287 * 0.00471 * -0.0000036732)
// newLat = Math.asin( 0.63456 + -0.000000013371 )
// newLat = Math.asin( 0.63455 )
// newLat = 0.68742

Now if we are calculating the North or South bearing – we can in fact stop there as we only need the latitude for that. While for East and West bearing we need only the longitude, you will have to calculate the latitude for these as it is required in the next step! The longitude is calculate by the original longitude in radians plus our Math.atan2 (as discussed in the previous blog post) result on sine of the bearing multiplied by sin of the distance multiplied by cosine of the original latitude, and cosine of the distance minus sine of the original latitude multiplied by sine of the new latitude. If you’re getting confused by this point – join the club! Let’s look at it in code as it makes more sense in my opinion!

// New longitude function lonR + atan2(sin(eRad) * sin(distance) * cos(latR), cos(distance) - sin(latR) * sin(newLat))

var newLon = lonR + Math.atan2(Math.sin(eRad) * Math.sin(distance) * Math.cos(latR), Math.cos(distance) - Math.sin(latR) * Math.sin(newLat));
// newLon = -1.97774 + Math.atan2(Math.sin(1.57080) * Math.sin(0.00471) * Math.cos(0.68745), Math.cos(0.00471) - Math.sin(0.68745) * Math.sin(0.68742))
// newLon = -1.9774 + Math.atan2(0.99999 * 0.004709 * 0.77287, 0.99998 - 0.63457 * 0.63455)
// newLon = -1.9774 + Math.atan2(0.00364, 0.99998 - 0.40267)
// Don't forget BODMAS!
// newLon = -1.9774 + Math.atan2(0.00364, 0.59731)
// Again as our X is greater than 0 atan2 because the inverse tangent of y over x rememberin that y is the first parameter
// newLon = -1.9774 + Math.atan(0.00609)
// newLon = -1.9774 + 0.00609
// newLon = -1.97131

Finally we convert our radians back to degrees by multiplying by 180 and dividing by pi

newLat = newLat * 180 / Math.PI;
// newLat = 39.3863
newLon = newLon * 180 / Math.PI;
// newLon = -112.94784

Right so let’s have a look on Google maps at 39.3863, -112.94784. Again I wouldn’t have posted this if it wasn’t close enough! Keeping in my the rounding (and I’ll admit sometimes just to show I haven’t forgotten to do an operation I round in a strange manner to show the number has changed!), I get 31.48 km due East of our starting point!

So how does this optimise getting our data? Well, I’m not going to work through the other 3 bearings needed – instead in my Google maps tab I’m just going to pull out rough values for North (39.660204,-113.3143567), South (39.118926, -113.318818) and West (39.393339, -113.666147). Again remember these aren’t exact figures – in production I would run the function without poor rounding for each of these four points. From here I essentially want to derive a square in my mind. I use the longitude (remember that’s counting from the Prime Meridian that runs through Greenwich, London, UK) on my West point, and the latitude (remember that’s counting from the equator) on my North point – these are my minimum values. I then take the longitude from my East point and latitude from my South point for my maximum values. If you imagine drawing a square from these minimum and maximum points (with the other corners being north latitude and east longitude and south latitude and west longitude), you have an area only slightly larger than the circle of our radius. Of course there are some people in the outer part of the arc – but we can always run our Haversine formula on the results from this square to finesse the results. For the most part, our results within this square will be correct, and because we now have numbers I can call the following function assuming use of MySQL

SELECT * FROM restaurants WHERE lat > minLatitude AND lat < maxLatitude AND lon > minLongitude AND lon < maxLongitude;

From this I receive only nearby results – if the requirement is strictly 30 kms I can then iterate over the results with the Haversine formula to ensure they are within the 30 kms radius – if the require is “nearby” and you are working on a smaller scale (say 5km), you could probably get away with just returning these results!

So there you have it! This sure caused me a headache to solve – hopefully this will help you with simplifying and optimising location related codes! Do you know of an easier way to do this? Start the conversation in the comments section over on our LinkedIn! We’d love to hear how you overcome location filtering!

function getPoint(startLat, startLon, bearingDeg, targetDistance) {
    var distance = targetDistance / 6371;
    var bearing = bearingDeg * Math.PI / 180;
    var latR = startLat * Math.PI / 180;
    var lonR = startLon * Math.PI / 180;
    var newLat = Math.asin(Math.sin(latR) * Math.cos(distance) + Math.cos(latR) * Math.sin(distance) * Math.cos(bearing));
    var newLon = lonR + Math.atan2(Math.sin(bearing) * Math.sin(distance) * Math.cos(latR), Math.cos(distance) - Math.sin(latR) * Math.sin(newLat));
    newLat = newLat * 180 / Math.PI;
    newLon = newLon * 180 / Math.PI;
    return newLat + ", " + newLon;
}