pg.ddx.io  pgsql-sql@postgresql.org mailing list archive  
help / color / mirror / Atom feed
Radius of a zip code
10+ messages / 6 participants
[nested] [flat]

* Radius of a zip code
@ 2003-12-26 23:42  Andy Lewis <jumboc@comcast.net>
  0 siblings, 2 replies; 10+ messages in thread

From: Andy Lewis @ 2003-12-26 23:42 UTC (permalink / raw)
  To: pgsql-sql

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

^ permalink  raw  reply  [nested|flat] 10+ messages in thread

* Re: Radius of a zip code
@ 2003-12-27 00:09  Tom Lane <tgl@sss.pgh.pa.us>
  parent: Andy Lewis <jumboc@comcast.net>
  1 sibling, 1 reply; 10+ messages in thread

From: Tom Lane @ 2003-12-27 00:09 UTC (permalink / raw)
  To: Andy Lewis <jumboc@comcast.net>; +Cc: pgsql-sql

"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



^ permalink  raw  reply  [nested|flat] 10+ messages in thread

* Re: Radius of a zip code
@ 2003-12-27 02:42  Michael Fuhr <mike@fuhr.org>
  parent: Andy Lewis <jumboc@comcast.net>
  1 sibling, 3 replies; 10+ messages in thread

From: Michael Fuhr @ 2003-12-27 02:42 UTC (permalink / raw)
  To: Andy Lewis <jumboc@comcast.net>; +Cc: pgsql-sql

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    |                 0
 75270   | Dallas        | TX    | 0.460576795779555
 75202   | Dallas        | TX    |  0.62326173788043
  .
  .
  .
 76012   | Arlington     | TX    |   19.644132573068
 75126   | Forney        | TX    |  19.8963253723536
 75024   | 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/



^ permalink  raw  reply  [nested|flat] 10+ messages in thread

* Re: Radius of a zip code
@ 2003-12-27 03:08  Joe Conway <mail@joeconway.com>
  parent: Michael Fuhr <mike@fuhr.org>
  2 siblings, 0 replies; 10+ messages in thread

From: Joe Conway @ 2003-12-27 03:08 UTC (permalink / raw)
  To: Michael Fuhr <mike@fuhr.org>; +Cc: Andy Lewis <jumboc@comcast.net>; pgsql-sql

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 ALIAS FOR $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 as z
     WHERE
      abs(z.lat - $lat1d) * 60 * 1.15078 <= $dist
      and abs(z.long - $lon1d) * 60 * 1.15078 <= $dist
      and zipdist($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





^ permalink  raw  reply  [nested|flat] 10+ messages in thread

* Re: Radius of a zip code
@ 2003-12-27 04:34  Bruno Wolff III <bruno@wolff.to>
  parent: Michael Fuhr <mike@fuhr.org>
  2 siblings, 1 reply; 10+ messages in thread

From: Bruno Wolff III @ 2003-12-27 04:34 UTC (permalink / raw)
  To: Michael Fuhr <mike@fuhr.org>; +Cc: Andy Lewis <jumboc@comcast.net>; pgsql-sql

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.



^ permalink  raw  reply  [nested|flat] 10+ messages in thread

* Re: Radius of a zip code
@ 2003-12-27 05:19  Michael Fuhr <mike@fuhr.org>
  parent: Bruno Wolff III <bruno@wolff.to>
  0 siblings, 1 reply; 10+ messages in thread

From: Michael Fuhr @ 2003-12-27 05:19 UTC (permalink / raw)
  To: pgsql-sql

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 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.

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

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



^ permalink  raw  reply  [nested|flat] 10+ messages in thread

* Re: Radius of a zip code
@ 2003-12-27 07:50  Bruno Wolff III <bruno@wolff.to>
  parent: Michael Fuhr <mike@fuhr.org>
  0 siblings, 1 reply; 10+ messages in thread

From: Bruno Wolff III @ 2003-12-27 07:50 UTC (permalink / raw)
  To: Michael Fuhr <mike@fuhr.org>; +Cc: pgsql-sql

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.



^ permalink  raw  reply  [nested|flat] 10+ messages in thread

* Re: Radius of a zip code
@ 2003-12-27 15:36  Andy Lewis <jumboc@comcast.net>
  parent: Michael Fuhr <mike@fuhr.org>
  2 siblings, 0 replies; 10+ messages in thread

From: Andy Lewis @ 2003-12-27 15:36 UTC (permalink / raw)
  To: 'Michael Fuhr' <mike@fuhr.org>; +Cc: pgsql-sql

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    |                 0
 75270   | Dallas        | TX    | 0.460576795779555
 75202   | Dallas        | TX    |  0.62326173788043
  .
  .
  .
 76012   | Arlington     | TX    |   19.644132573068
 75126   | Forney        | TX    |  19.8963253723536
 75024   | 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/




^ permalink  raw  reply  [nested|flat] 10+ messages in thread

* Re: Radius of a zip code
@ 2003-12-27 18:21  Tom Lane <tgl@sss.pgh.pa.us>
  parent: Bruno Wolff III <bruno@wolff.to>
  0 siblings, 0 replies; 10+ messages in thread

From: Tom Lane @ 2003-12-27 18:21 UTC (permalink / raw)
  To: Bruno Wolff III <bruno@wolff.to>; +Cc: Michael Fuhr <mike@fuhr.org>; pgsql-sql

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



^ permalink  raw  reply  [nested|flat] 10+ messages in thread

* Re: Radius of a zip code
@ 2004-01-01 04:56  Randolf Richardson <rr@8x.ca>
  parent: Tom Lane <tgl@sss.pgh.pa.us>
  0 siblings, 0 replies; 10+ messages in thread

From: Randolf Richardson @ 2004-01-01 04:56 UTC (permalink / raw)
  To: pgsql-sql

"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.



^ permalink  raw  reply  [nested|flat] 10+ messages in thread


end of thread, other threads:[~2004-01-01 04:56 UTC | newest]

Thread overview: 10+ messages (download: mbox mbox.gz follow: Atom feed)
-- links below jump to the message on this page --
2003-12-26 23:42 Radius of a zip code Andy Lewis <jumboc@comcast.net>
2003-12-27 00:09 ` Tom Lane <tgl@sss.pgh.pa.us>
2004-01-01 04:56   ` Randolf Richardson <rr@8x.ca>
2003-12-27 02:42 ` Michael Fuhr <mike@fuhr.org>
2003-12-27 03:08   ` Joe Conway <mail@joeconway.com>
2003-12-27 04:34   ` Bruno Wolff III <bruno@wolff.to>
2003-12-27 05:19     ` Michael Fuhr <mike@fuhr.org>
2003-12-27 07:50       ` Bruno Wolff III <bruno@wolff.to>
2003-12-27 18:21         ` Tom Lane <tgl@sss.pgh.pa.us>
2003-12-27 15:36   ` Andy Lewis <jumboc@comcast.net>

This inbox is served by DDX for PostgreSQL; see mirroring instructions
for how to clone and mirror all data and code used for this inbox