X-Original-To: pgsql-sql-postgresql.org@localhost.postgresql.org
Received: from localhost (neptune.hub.org [200.46.204.2])
	by svr1.postgresql.org (Postfix) with ESMTP id 91C09D1B471
	for <pgsql-sql-postgresql.org@localhost.postgresql.org>;
	Sat, 27 Dec 2003 03:13:29 +0000 (GMT)
Received: from svr1.postgresql.org ([200.46.204.71])
	by localhost (neptune.hub.org [200.46.204.2]) (amavisd-new, port 10024)
	with ESMTP id 70233-10
	for <pgsql-sql-postgresql.org@localhost.postgresql.org>;
	Fri, 26 Dec 2003 23:12:44 -0400 (AST)
Received: from joeconway.com (66-146-172-86.skyriver.net [66.146.172.86])
	by svr1.postgresql.org (Postfix) with ESMTP id E8D37D1B448
	for <pgsql-sql@postgresql.org>; Fri, 26 Dec 2003 23:12:39 -0400 (AST)
Received: from [192.168.5.3] (account jconway HELO joeconway.com)
	by joeconway.com (CommuniGate Pro SMTP 4.1.8)
	with ESMTP-TLS id 746616; Fri, 26 Dec 2003 19:09:40 -0800
Message-ID: <3FECF7A3.2050108@joeconway.com>
Date: Fri, 26 Dec 2003 19:08:19 -0800
From: Joe Conway <mail@joeconway.com>
User-Agent: Mozilla/5.0 (Windows; U; Windows NT 5.0; en-US;
	rv:1.5) Gecko/20031007
X-Accept-Language: en-us, en
MIME-Version: 1.0
To: Michael Fuhr <mike@fuhr.org>
Cc: Andy Lewis <jumboc@comcast.net>, pgsql-sql@postgresql.org
Subject: Re: Radius of a zip code
References: <000b01c3cc09$e043fd70$0201a8c0@andy2>
	<20031226194244.A15916@quality.qadas.com>
In-Reply-To: <20031226194244.A15916@quality.qadas.com>
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit
X-Virus-Scanned: by amavisd-new at postgresql.org
X-Archive-Number: 200312/193
X-Sequence-Number: 16359

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