I was working on an obscure but very interesting coding problem this week. It turned out to have an interesting solution. The problem is a bit tricky to explain. I was working with latitude-longitude data but the idea works with any real-valued variables. I’ll illustrate with the latitude part. A latitude is a value between -90.0 and +90.0 that is the “up-down” location on a map. Latitude +90.0 is the North Pole and latitude 0 is the equator.

For reasons that are a very long story (see http://msdn.microsoft.com/en-us/magazine/jj133823.aspx) I had a set of 1,800 latitude intervals of width 0.1 degree (roughly 7 miles or 11 kilometers) and each interval has an index/ID:

Interval Index ---------------------- [-90.0, -89.9) 0 [-89.9, -89.8) 1 [-89.8, -89.7) 2 . . . [ 0.0, +0.1) 900 [ +0.1, +0.2) 901 . . . [ 89.7, 89.8) 1797 [ 89.8, 89.9) 1798 [ 89.9, 90.0] 1799

Notice the left part of each interval has a square bracket, meaning greater than or equal to, but the right part of each interval, except for the last interval, has a parenthesis, meaning strictly less than.

The problem I needed to solve was to write a function that accepts a latitude value such as 89.7321 and returns the correct index, 1797.

The first part of the solution is to map the input latitude to the left part its interval. The second part is to map the left part of the interval to the index.

Here is C# code that maps a latitude, lat, to the left part of the interval:

double lat = 89.7321; // left interval is 89.7 double timesTen = lat * 10; // 897.321 double floorOfTimesTen = Math.Floor(timesTen); // 897.0 double intervalLeft = floorOfTimesTen / 10.0; // 89.7

Notice I use the Floor function rather than casting to an int because latitudes can be negative. The second part essentially solves the equation of a line, y = mx + b, where y is the desired index and x is the input latitude. If the first interval-index pair is point p1 then (x1,y1) = (-90.0, 0). And then p2 = (x2,y2) = (+89.9, 1799). Solving for the slope: m = (y2-y1) / (x2-1) = (1799 – 0) / 89.9 – (-90.0) = 10. Now using p1, we can get the y-intercept: b = y – mx = 0 – (10)(-90.0) = 900.

So the second part of the code looks like this:

return (int)(10 * intervalLeft) + 900;

Putting this all together, and dealing with the last interval special case, we get:

static int LatIndex(double lat) { if (lat == 90.0) return 1799; // special case double floorOfTimesTen = Math.Floor(lat * 10); double intervalLeft = floorOfTimesTen / 10.0; return (int)(10 * intervalLeft) + 900; }

Very interesting problem. I use this function, and its cousin LonIndex, to compute a sector ID for a latitude-longitude. I have sector IDs stored in a database so that when I want to find geo-location data points that are near to each other, I just look up data points that have the same sector ID, which, for my large database, is literally hundreds of thousands of times faster than computing distances.