Haversine Formula unit errors - PL/SQL

Question

Overall goal: Find all zip codes within 30 miles of a given zip.

Data Source: List of Zip Codes, with Latitude and Longitude as 'x' & 'y'.

Example:

create or replace view zips as
select 37171 zip, 36.362704 y, -87.30434 x,'Southside' City, 'TN' State from dual 
union
select 37212, 36.133012, -86.802764, 'Nashville', 'TN' from dual 
union
select 37027, 36.00245, -86.791159, 'Brentwood', 'TN' from dual 
union
select 37191, 36.501196, -87.539057, 'Woodlawn', 'TN' from dual 
union
select 37067, 35.928406, -86.805538, 'Franklin', 'TN' from dual ;

What I've tried: I found something that looked like it was going to solve all my problems:

CREATE OR REPLACE FUNCTION distance (Lat1 IN NUMBER,
                                     Lon1 IN NUMBER,
                                     Lat2 IN NUMBER,
                                     Lon2 IN NUMBER,
                                     Radius IN NUMBER DEFAULT 3963) RETURN NUMBER IS
 -- Convert degrees to radians
 DegToRad NUMBER := 57.29577951;

BEGIN
  RETURN(NVL(Radius,0) * ACOS((sin(NVL(Lat1,0) / DegToRad) * SIN(NVL(Lat2,0) / DegToRad)) +
        (COS(NVL(Lat1,0) / DegToRad) * COS(NVL(Lat2,0) / DegToRad) *
         COS(NVL(Lon2,0) / DegToRad - NVL(Lon1,0)/ DegToRad))));
END;

However, it was giving me funky numbers - consistently too small (I know it's not 1 mile from Franklin to Nashville), but not by a constant factor that I could determine. (Code I used to test is below - I just compared distances on Google Maps with these distances)

select b.zip
  ,b.city
  ,b.state
  ,distance(a.x,a.y,b.x,b.y) distance
from zips a, zips b
where a.zip=37067
order by distance;

So, I thought maybe the geodata group had a different way of recording lat & long and I found Wikipedia's Haversine formula and turned that into a function

CREATE OR REPLACE FUNCTION distance (Lat1_d IN NUMBER,
                                     Lon1_d IN NUMBER,
                                     Lat2_d IN NUMBER,
                                     Lon2_d IN NUMBER,
                                     Radius IN NUMBER DEFAULT 3959) RETURN NUMBER IS
 -- Convert degrees to radians
 DegToRad NUMBER := .0174532925;
 Lat1 NUMBER := Lat1_d * DegToRad;
 Lon1 NUMBER := Lon1_d * DegToRad;
 Lat2 NUMBER := Lat2_d * DegToRad;
 Lon2 NUMBER := Lon2_d * DegToRad;

BEGIN

  RETURN 2*Radius * 
    asin(
      sqrt(
        power(sin((Lat2-Lat1)/2),2) + cos(Lat1) * cos(Lat2) * power(sin((Lon2-Lon1)/2),2)
      )
    )
        ;
END;

Same issue - results are nearly exactly the same, which leaves me wondering what the problem is. What transformation am I not doing? What am I missing?

Note that I would be glad to assume a flat earth and use Pythagorean distances since I'm only going to look at 30-50 miles for each set of zip codes I look at, but I'd need to figure out what about my latitude / longitude conversion is going wrong there as well.

Answer 1

Looks like a parameter issue. The header is:

CREATE OR REPLACE FUNCTION distance (Lat1 IN NUMBER,
                                     Lon1 IN NUMBER,
                                     Lat2 IN NUMBER,
                                     Lon2 IN NUMBER,

And you're calling it as follows:

distance(a.x,a.y,b.x,b.y)

But your view is:

create or replace view zips as
select 37171 zip, 36.362704 y, -87.30434 x,'Southside' City, 'TN' State from dual 
union
select 37212, 36.133012, -86.802764, 'Nashville', 'TN' from dual 
union
select 37027, 36.00245, -86.791159, 'Brentwood', 'TN' from dual 
union
select 37191, 36.501196, -87.539057, 'Woodlawn', 'TN' from dual 
union
select 37067, 35.928406, -86.805538, 'Franklin', 'TN' from dual ;

According to the view, x is the longitude and y is the latitude. So the function call should be:

distance(a.y,a.x,b.y,b.x)

Alternatively, you can swap the "x" and "y" aliases in the view and leave everything else as is.