I have created a table where I store Lat/Long. data along with ZIP codes. Here is my procedure to find ZIPs within certain radius:
CREATE PROCEDURE geodistZCTA (IN userid char(5), IN dist double)
BEGIN
declare mylon double;
declare mylat double;
declare lon1 float;
declare lon2 float;
declare lat1 float;
declare lat2 float;
-- get the original lon and lat for the userid
select Y(location), X(location) into mylon, mylat from zcta where zip=userid limit 1;
-- calculate lon and lat for the rectangle:
set lon1 = mylon-dist/abs(cos(radians(mylat))*69);
set lon2 = mylon+dist/abs(cos(radians(mylat))*69);
set lat1 = mylat-(dist/69);
set lat2 = mylat+(dist/69);
-- run the query:
-- select astext(location), zip from zcta where st_within(location, envelope(linestring(point(lon1, lat1), point(lon2, lat2)))) order by st_distance(point(mylon, mylat), location) limit 10;
-- select zip,(3956 * 2 * ASIN(SQRT(POWER(SIN((lat1 - abs(lat2)) * pi()/180 / 2), 2) + COS(abs(lat1) * pi()/180 ) * COS(abs(lat2) * pi()/180) * POWER(SIN((lon1 - lon2) * pi()/180 / 2), 2) ))) AS distance FROM zcta HAVING distance<=dist ORDER BY distance ASC;
-- SELECT zip,X(location),Y(location),(((ACOS(SIN(mylat * PI() / 180) * SIN(X(location) * PI() / 180) + COS(mylat * PI() / 180) * COS(X(location) * PI() / 180) * COS((Y(location)-mylon) * PI() / 180)) * 180 / PI()) * 60 * 1.1515)) AS distance FROM zcta HAVING distance<=dist ORDER BY distance ASC;
select zip,X(location),Y(location),
( 3959 * acos( cos( radians(mylat) )
* cos( radians( X(location) ) )
* cos( radians( Y(location) ) - radians(mylon) )
+ sin( radians(mylat) )
* sin( radians( X(location) ) ) ) ) AS distance
FROM zcta
where X(location) between lat1 and lat2
and Y(location) between lon1 and lon2
HAVING distance <= dist ORDER BY distance ASC;
END;
Now when I call the procedure like:
call geodistZCTA('03804',5);
I get the following result:
# zip, X(location), Y(location), distance
-------------------------------------------------------
'03042', '43.045076', '-71.07095', '3.9798046354127266'
But when I measure the distance using http://www.allplaces.us/dbz.cgi I get 16.2 miles ( 26.1 km )
Any idea of what I am doing wrong here?