Re: Radius of a zip code

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