Re: Radius of a zip code

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

"Andy Lewis" <jumboc(at)comcast(dot)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

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/

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

On Fri, Dec 26, 2003 at 19:42:44 -0700,
Michael Fuhr <mike(at)fuhr(dot)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.

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(at)fuhr(dot)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/

On Fri, Dec 26, 2003 at 22:19:52 -0700,
Michael Fuhr <mike(at)fuhr(dot)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.

Thanks All for your suggestions, I have enough information to construct
what I need.

-----Original Message-----
From: Michael Fuhr [mailto:mike(at)fuhr(dot)org]
Sent: Friday, December 26, 2003 8:43 PM
To: Andy Lewis
Cc: pgsql-sql(at)postgresql(dot)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/

Bruno Wolff III <bruno(at)wolff(dot)to> writes:
> Michael Fuhr <mike(at)fuhr(dot)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

"tgl(at)sss(dot)pgh(dot)pa(dot)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(at)8x(dot)ca
Vancouver, British Columbia, Canada

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