Geolocation SQL query not finding exact location

Question

I have been testing my geolocation query for some time now and I haven't found any issues with it until now.

I am trying to search for all cities within a given radius, often times I'm searching for cities surrounding a city using that city's coords, but recently I tried searching around a city and found that the city itself was not returned.

I have these cities as an excerpt in my database:

city            latitude    longitude  
Saint-Mathieu   45.316708   -73.516253  
Saint-Édouard   45.233374   -73.516254  
Saint-Michel    45.233374   -73.566256  
Saint-Rémi      45.266708   -73.616257  

But when I run my query around the city of Saint-Rémi, with the following query...

SELECT tblcity.city, tblcity.latitude, tblcity.longitude, 
truncate((degrees(acos( sin(radians(tblcity.latitude)) 
* sin(radians(45.266708)) 
+ cos(radians(tblcity.latitude)) 
* cos(radians(45.266708)) 
* cos(radians(tblcity.longitude - -73.616257) ) ) ) 
* 69.09*1.6),1) as distance 
FROM tblcity HAVING distance < 10 ORDER BY distance desc 

I get these results:

city            latitude    longitude     distance  
Saint-Mathieu   45.316708   -73.516253    9.5  
Saint-Édouard   45.233374   -73.516254    8.6  
Saint-Michel    45.233374   -73.566256    5.3  

The town of Saint-Rémi is missing from the search.

So I tried a modified query hoping to get a better result:

SELECT tblcity.city, tblcity.latitude, tblcity.longitude, 
truncate(( 6371 * acos( cos( radians( 45.266708 ) ) 
* cos( radians( tblcity.latitude ) ) 
* cos( radians( tblcity.longitude ) 
- radians( -73.616257 ) ) 
+ sin( radians( 45.266708 ) ) 
* sin( radians( tblcity.latitude ) ) ) ),1) AS distance 
FROM tblcity HAVING distance < 10 ORDER BY distance desc 

But I get the same result...

However, if I modify Saint-Rémi's coords slighly by changing the last digit of the lat or long by 1, both queries will return Saint-Rémi. Also, if I center the query on any of the other cities above, the searched city is returned in the results.

Can anyone shed some light on what may be causing my queries above to not display the searched city of Saint-Rémi? I have added a sample of the table (with extra fields removed) below.

I'm using MySQL 5.0.45, thanks in advance.

CREATE TABLE `tblcity` ( 
`IDCity` int(1) NOT NULL auto_increment, 
`City` varchar(155) NOT NULL default '', 
`Latitude` decimal(9,6) NOT NULL default '0.000000', 
`Longitude` decimal(9,6) NOT NULL default '0.000000', 
PRIMARY KEY (`IDCity`) 
) ENGINE=MyISAM AUTO_INCREMENT=52743 DEFAULT CHARSET=latin1 AUTO_INCREMENT=52743; 

INSERT INTO `tblcity` (`city`, `latitude`, `longitude`) VALUES 
('Saint-Mathieu', 45.316708, -73.516253), 
('Saint-Édouard', 45.233374, -73.516254), 
('Saint-Michel', 45.233374, -73.566256), 
('Saint-Rémi', 45.266708, -73.616257); 

Answer 1

In your first query, I believe you've inverted the longitudes in the subtraction. The Spherical Law of Cosines is:

d = acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(long2−long1))*R

If lat1 is substituted with tblcity.latitude, long1 must be substituted with tblcity.longitude. I think you've accidentally substituted long2 in your query. Does this one work better?

SELECT tblcity.city, tblcity.latitude, tblcity.longitude, 
truncate((degrees(acos( sin(radians(tblcity.latitude)) 
* sin(radians(45.266708)) 
+ cos(radians(tblcity.latitude)) 
* cos(radians(45.266708)) 
* cos(radians(-73.616257 - tblcity.longitude) ) ) ) 
* 69.09*1.6),1) as distance 
FROM tblcity HAVING distance < 10 ORDER BY distance desc 

I haven't looked into your second query yet, but hopefully that helps.

Answer 2

You are using the "spherical law of cosines" formula, which is susceptable to rounding error at small distances (like zero!) -- see this discussion. The long expression that is fed into acos() is evaluating to slightly more than 1.0, which is out of bounds.

Here's the problem illustrated using Python to do the calculations:

>>> from math import sin, cos, acos, radians
>>> lat = radians(45.266708)
>>> long_expression = sin(lat) * sin(lat) + cos(lat) * cos(lat) * cos(0.0)
>>> repr(long_expression)
'1.0000000000000002'
>>> acos(long_expression)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: math domain error
>>>

It seems that MySQL is substituting NULL instead of raising an exception. I know little about MySQL, but you should be able to overcome that by doing something like ifnull(acos(long_expression), 0.0) or coalesce(acos(long_expression), 0.0).

Alternatively you could use the haversine formula, which shifts the rounding problem from your door-step to the opposite side of the earth.

Update: I've tested using that formula in Python to calculate the should-be-zero distance between a point and the same point, for each of the 37582 unique (lat, lon) 2-tuples in a file of US zip codes.

Of these:

• 31591 (84.1%) produced a zero distance
• 4244 (11.3%) produced a distance of 9.5 cm.
• 831 (2.2%) produced a distance of 13.4 cm.
• 916 (2.4%) produced a "cos" value of 1.0000000000000002 which would cause an exception in acos() if not detected and avoided.

It appears that explicit testing for lat1 = lat2 and lon1 = lon2 and avoiding the formula in that case (just use zero) might be a good idea -- it would give a consistent answer and avoid puzzlement.