Postgres Pro
Downloads
Home   >   mailing lists

Thread: Radius of a zip code

Radius of a zip code

From
"Andy Lewis"
Date:
Hello all,
 
I was trying to find all zip codes within a given zip code or radius.
 
I have map points and Latitude and Longitude in my zip table.
 
I remember seeing a post or two referencing this but can't see to find it.
 
I've tried the following with no luck:
 
-- 20 Miles
--select 20 * 360.0 / (7900.0 * atan2(1,1) * 4.0);
select * from zip_code where map_loc @ circle(map_point('dallas','tx','75201'), .290105212724467 ) order by city
 
Anyone that has this experience, can you validate this for correctness?
 
Thanks in advance,
 
Andy

Re: Radius of a zip code

From
Tom Lane
Date:
"Andy Lewis" <jumboc@comcast.net> writes:
> I was trying to find all zip codes within a given zip code or radius.

I think there are canned solutions for this available in PostGIS ---
have you looked at that?

> I've tried the following with no luck:
> -- 20 Miles
> --select 20 * 360.0 / (7900.0 * atan2(1,1) * 4.0);
> select * from zip_code where map_loc @
> circle(map_point('dallas','tx','75201'), .290105212724467 ) order by
> city

I'm guessing that the big problem is that you didn't measure longitude
and latitude in identical units in your table, so your "circle" isn't
real circular, and the smaller problem is that "miles" converts to
"degrees of arc" differently at different latitudes.
        regards, tom lane

Re: Radius of a zip code

From
Michael Fuhr
Date:
On Fri, Dec 26, 2003 at 05:42:08PM -0600, Andy Lewis wrote:
> I was trying to find all zip codes within a given zip code or radius.
>  
> I have map points and Latitude and Longitude in my zip table.
>  
> I remember seeing a post or two referencing this but can't see to find
> it.

The code in contrib/earthdistance in the PostgreSQL source code might
be what you're looking for.  I haven't used it myself, as I had already
written a function I needed for another DBMS and ported it to PostgreSQL.

> I've tried the following with no luck:
>  
> -- 20 Miles
> --select 20 * 360.0 / (7900.0 * atan2(1,1) * 4.0);
> select * from zip_code where map_loc @
> circle(map_point('dallas','tx','75201'), .290105212724467 ) order by
> city

This isn't related to the problem, but is there a reason your map_point
function requires city, state, and zip code?  If you know the zip code
then you shouldn't need the city and state.

> Anyone that has this experience, can you validate this for correctness?

I have several databases with lat/lon coordinates and frequently make
"show me all records within a certain distance of this point" queries.
I wrote a haversine() function that uses the Haversine Formula to
calculate the great circle distance between two points on a sphere
(assuming the earth is a perfect sphere is accurate enough for my uses).
Here's a web site with related info:

http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1

Here's an example of how I use the haversine() function.  I'm not using
PostgreSQL's geometric types -- latitude and longitude are stored in
separate fields.  The function takes two lat/lon coordinates in degrees
and optionally a radius (the default is 3956.0, the approximate radius
of the earth in miles); it returns the distance in whatever units the
radius is in.

SELECT a.zipcode, a.city, a.state,      haversine(a.latitude, a.longitude, b.latitude, b.longitude) AS dist
FROM zipcode AS a, zipcode AS b
WHERE b.zipcode = 75201 AND haversine(a.latitude, a.longitude, b.latitude, b.longitude) <= 20
ORDER BY dist;
zipcode |     city      | state |       dist        
---------+---------------+-------+-------------------75201   | Dallas        | TX    |                 075270   |
Dallas       | TX    | 0.46057679577955575202   | Dallas        | TX    |  0.62326173788043 . . .76012   | Arlington
| TX    |   19.64413257306875126   | Forney        | TX    |  19.896325372353675024   | Plano         | TX    |
19.9884653971924
(106 rows)

As for validating the function's correctness, I'm using a well-known
formula and I've compared the function's output to distances measured on
a map.  I wouldn't use it for missile targeting, but it's sufficiently
accurate for "show me all stores within 20 miles of my home."

Here's the meat of the function (written in C); the coordinates have by
now been converted to radians:
 dlat = lat2 - lat1; dlon = lon2 - lon1;
 a1 = sin(dlat / 2.0); a2 = sin(dlon / 2.0);
 a = (a1 * a1) + cos(lat1) * cos(lat2) * (a2 * a2); c = 2.0 * atan2(sqrt(a), sqrt(1.0 - a));
 dist = radius * c;

If anybody's interested I'll post the entire file.

-- 
Michael Fuhr
http://www.fuhr.org/~mfuhr/

Re: Radius of a zip code

From
Joe Conway
Date:
Michael Fuhr wrote:
> I wrote a haversine() function that uses the Haversine Formula to
> calculate the great circle distance between two points on a sphere
> (assuming the earth is a perfect sphere is accurate enough for my uses).
> Here's a web site with related info:
> 
> http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1

[...snip...]

> Here's the meat of the function (written in C); the coordinates have by
> now been converted to radians:

[...snip...]

> If anybody's interested I'll post the entire file.

FWIW, here's a plpgsql function I wrote a while ago based on the 
Haversine formula:

CREATE FUNCTION "zipdist" (float8,float8,float8,float8 ) RETURNS float8 AS '   DECLARE     lat1 ALIAS FOR $1;     lon1
ALIASFOR $2;     lat2 ALIAS FOR $3;     lon2 ALIAS FOR $4;     dist float8;   BEGIN     dist := 0.621 * 6371.2 * 2 *
  atan2( sqrt(abs(0 + pow(sin(radians(lat2)/2 -       radians(lat1)/2),2) + cos(radians(lat1)) * cos(radians(lat2)) *
   pow(sin(radians(lon2)/2 - radians(lon1)/2),2))),sqrt(abs(1 -       pow(sin(radians(lat2)/2 - radians(lat1)/2),2) +
cos(radians(lat1))*       cos(radians(lat2)) * pow(sin(radians(lon2)/2 - 
 
radians(lon1)/2),2))));     return dist;   END;
' LANGUAGE 'plpgsql';

I used the following PHP code to start looking for a match in a small 
circle, and then expand it if no matches were found:

$dist = INIT_DIST;
$cnt = 0;
$cntr = 0;
do {     if ((! $zip == "") && (! $dist <= 0)) {         $sql = get_zip_sql($lon1d,$lat1d,$dist,$numtoshow);
$rs= connexec($conn,$sql);         $rsf = rsfetchrs($rs);         $dist *= 2;         $cntr++;     } else {
$cntr= 10;     }
 
} while (count($rsf) < $numadvisorstoshow && $cntr < 10);

Hopefully you get the idea. You can narrow the results using a box to 
make the query perform better, and then sort by distance to get the 
closest alternative. Here's the related part of get_zip_sql():

function get_zip_sql($lon1d,$lat1d,$dist,$numtoshow)
{  $sql = "    SELECT DISTINCT <fields>    FROM tbl_a AS a     ,tbl_d AS d     ,tbl_a_zipcodes AS az     ,tbl_zipcodes
asz    WHERE     abs(z.lat - $lat1d) * 60 * 1.15078 <= $dist     and abs(z.long - $lon1d) * 60 * 1.15078 <= $dist
andzipdist($lat1d,$lon1d,lat,long) <= $dist     and z.zip = az.zipcode     <other criteria>    ORDER BY    LIMIT
$numtoshow; ";
 
  return $sql;
}

The "X * 60 * 1.15078" converts differences in degrees lat/long into 
rough distances in miles.

Hope this helps.

Joe

Re: Radius of a zip code

From
Bruno Wolff III
Date:
On Fri, Dec 26, 2003 at 19:42:44 -0700, Michael Fuhr <mike@fuhr.org> wrote:
> 
> I have several databases with lat/lon coordinates and frequently make
> "show me all records within a certain distance of this point" queries.
> I wrote a haversine() function that uses the Haversine Formula to
> calculate the great circle distance between two points on a sphere
> (assuming the earth is a perfect sphere is accurate enough for my uses).
> Here's a web site with related info:

The distance operator in contrib/earthdistance got changed to use
haversine instead of the naive formula in 7.3. In 7.4 it also provides
some functions that work with contrib/cube that allow for indexed
searches.

Re: Radius of a zip code

From
Michael Fuhr
Date:
On Fri, Dec 26, 2003 at 10:34:04PM -0600, Bruno Wolff III wrote:
> On Fri, Dec 26, 2003 at 19:42:44 -0700,
>   Michael Fuhr <mike@fuhr.org> wrote:
> > 
> > I have several databases with lat/lon coordinates and frequently make
> > "show me all records within a certain distance of this point" queries.
> > I wrote a haversine() function that uses the Haversine Formula to
> > calculate the great circle distance between two points on a sphere
> > (assuming the earth is a perfect sphere is accurate enough for my uses).
> > Here's a web site with related info:
> 
> The distance operator in contrib/earthdistance got changed to use
> haversine instead of the naive formula in 7.3. In 7.4 it also provides
> some functions that work with contrib/cube that allow for indexed
> searches.

I'll have to take a closer look at contrib/earthdistance.  I'm using the
function I wrote for legacy reasons -- I had ported an application from
another DBMS to PostgreSQL and wanted to make as few changes as possible,
so I ported the haversine() function that I had already written.

Incidentally, I see the following in README.earthdistance:
 A note on testing C extensions - it seems not enough to drop a function and re-create it - if I change a function, I
haveto stop and restart the backend for the new version to be seen. I guess it would be too messy to track which
functionsare added from a .so and do a dlclose when the last one is dropped.
 

Maybe you've already figured it out, but LOAD should allow you to reload
a .so file without having to restart the backend.

http://www.postgresql.org/docs/current/static/sql-load.html

-- 
Michael Fuhr
http://www.fuhr.org/~mfuhr/

Re: Radius of a zip code

From
Bruno Wolff III
Date:
On Fri, Dec 26, 2003 at 22:19:52 -0700, Michael Fuhr <mike@fuhr.org> wrote:
> 
> Incidentally, I see the following in README.earthdistance:
> 
>   A note on testing C extensions - it seems not enough to drop a function
>   and re-create it - if I change a function, I have to stop and restart
>   the backend for the new version to be seen. I guess it would be too
>   messy to track which functions are added from a .so and do a dlclose
>   when the last one is dropped.
> 
> Maybe you've already figured it out, but LOAD should allow you to reload
> a .so file without having to restart the backend.


I didn't write that. It came from the person(s) who worked on earthdistance
before I touched it.

Re: Radius of a zip code

From
"Andy Lewis"
Date:
Thanks All for your suggestions, I have enough information to construct
what I need.

-----Original Message-----
From: Michael Fuhr [mailto:mike@fuhr.org] 
Sent: Friday, December 26, 2003 8:43 PM
To: Andy Lewis
Cc: pgsql-sql@postgresql.org
Subject: Re: [SQL] Radius of a zip code


On Fri, Dec 26, 2003 at 05:42:08PM -0600, Andy Lewis wrote:
> I was trying to find all zip codes within a given zip code or radius.
>  
> I have map points and Latitude and Longitude in my zip table.
>  
> I remember seeing a post or two referencing this but can't see to find

> it.

The code in contrib/earthdistance in the PostgreSQL source code might be
what you're looking for.  I haven't used it myself, as I had already
written a function I needed for another DBMS and ported it to
PostgreSQL.

> I've tried the following with no luck:
>  
> -- 20 Miles
> --select 20 * 360.0 / (7900.0 * atan2(1,1) * 4.0);
> select * from zip_code where map_loc @ 
> circle(map_point('dallas','tx','75201'), .290105212724467 ) order by 
> city

This isn't related to the problem, but is there a reason your map_point
function requires city, state, and zip code?  If you know the zip code
then you shouldn't need the city and state.

> Anyone that has this experience, can you validate this for 
> correctness?

I have several databases with lat/lon coordinates and frequently make
"show me all records within a certain distance of this point" queries. I
wrote a haversine() function that uses the Haversine Formula to
calculate the great circle distance between two points on a sphere
(assuming the earth is a perfect sphere is accurate enough for my uses).
Here's a web site with related info:

http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1

Here's an example of how I use the haversine() function.  I'm not using
PostgreSQL's geometric types -- latitude and longitude are stored in
separate fields.  The function takes two lat/lon coordinates in degrees
and optionally a radius (the default is 3956.0, the approximate radius
of the earth in miles); it returns the distance in whatever units the
radius is in.

SELECT a.zipcode, a.city, a.state,      haversine(a.latitude, a.longitude, b.latitude, b.longitude) AS
dist FROM zipcode AS a, zipcode AS b WHERE b.zipcode = 75201 AND haversine(a.latitude, a.longitude, b.latitude,
b.longitude)<= 20
 
ORDER BY dist;
zipcode |     city      | state |       dist        
---------+---------------+-------+-------------------75201   | Dallas        | TX    |                 075270   |
Dallas       | TX    | 0.46057679577955575202   | Dallas        | TX    |  0.62326173788043 . . .76012   | Arlington
| TX    |   19.64413257306875126   | Forney        | TX    |  19.896325372353675024   | Plano         | TX    |
19.9884653971924
(106 rows)

As for validating the function's correctness, I'm using a well-known
formula and I've compared the function's output to distances measured on
a map.  I wouldn't use it for missile targeting, but it's sufficiently
accurate for "show me all stores within 20 miles of my home."

Here's the meat of the function (written in C); the coordinates have by
now been converted to radians:
 dlat = lat2 - lat1; dlon = lon2 - lon1;
 a1 = sin(dlat / 2.0); a2 = sin(dlon / 2.0);
 a = (a1 * a1) + cos(lat1) * cos(lat2) * (a2 * a2); c = 2.0 * atan2(sqrt(a), sqrt(1.0 - a));
 dist = radius * c;

If anybody's interested I'll post the entire file.

-- 
Michael Fuhr
http://www.fuhr.org/~mfuhr/

Re: Radius of a zip code

From
Tom Lane
Date:
Bruno Wolff III <bruno@wolff.to> writes:
>   Michael Fuhr <mike@fuhr.org> wrote:
>> Maybe you've already figured it out, but LOAD should allow you to reload
>> a .so file without having to restart the backend.

> I didn't write that. It came from the person(s) who worked on earthdistance
> before I touched it.

I've removed the incorrect comment.
        regards, tom lane

Re: Radius of a zip code

From
Randolf Richardson
Date:
"tgl@sss.pgh.pa.us (Tom Lane)" wrote in comp.databases.postgresql.sql:

[sNip]
> I'm guessing that the big problem is that you didn't measure longitude
> and latitude in identical units in your table, so your "circle" isn't
> real circular, and the smaller problem is that "miles" converts to
> "degrees of arc" differently at different latitudes.
       Don't forget that there are two different types of "miles" which need 
to be considered when measuring distances:
               1 statute/land mile = 1.609 km               1 nautical/sea mile = 1.85 km
       Since kilometers are consistent over land and water (and in the great 
vacuum of space), the metric system should always be used to ensure clarity, 
unless the only land masses the user is concerned with have no bodies of 
water.
       =)

-- 
Sir Randolf, noble spam fighter - rr@8x.ca
Vancouver, British Columbia, Canada

Please do not eMail me directly when responding
to my postings in the newsgroups.