[postgis-users] -- Still having trouble retrieving properties 100 meters around a point --

[Michael Fuhr]

On Tue, Sep 19, 2006 at 09:23:52PM -0700, Rafael Orta wrote:
> create table smarter.vow_properties (property_id serial, address varchar(150), latitude numeric (7,4) , longitude numeric (7,4))

As Greg pointed out, unless you need the lat/lon columns then you
could omit them and use only a geometry column.

CREATE TABLE smarter.vow_properties (
    property_id  serial PRIMARY KEY,
    address      text
);

> select AddGeometryColumn('smarter','vow_properties','geom',32011,'POINT',2);

SRID 32011 is NAD27 and feet; if you want meters then you could use
32618 (WGS84 / UTM zone 18N, which covers New Jersey).

SELECT AddGeometryColumn('smarter', 'vow_properties', 'geom', 32618, 'POINT', 2);

> insert into smarter.vow_properties(address, latitude, longitude) values ('12 LONG BRIDGE DR,MOUNT LAUREL, NJ 08054','39.940252','-074.921916');

Rather than inserting rows and then updating them, you could create
the geometries during the insert.  If we assume that the lat/lon
coordinates are WGS84 (SRID 4326) and that you want meters (SRID
32618) and that you don't need the lat/lon columns, then you could
do the inserts like this:

INSERT INTO smarter.vow_properties (address, geom)
  VALUES ('12 LONG BRIDGE DR, MOUNT LAUREL, NJ 08054',
	  Transform(SetSRID(MakePoint(-74.921916, 39.940252), 4326), 32618));

Another way to write the Transform would be:

Transform(GeomFromEWKT('SRID=4326;POINT(-74.921916 39.940252)'), 32618)

After the inserts, create an index on the geometry column and update
the query planner's statistics (you could skip creating the index
if the table doesn't have many rows).

CREATE INDEX vow_properties_geom_idx ON smarter.vow_properties USING gist (geom);
ANALYZE smarter.vow_properties;

> The search below suppose to retrieve properties 100 meters around the point
> logitude = -074.905868 latitude = 39.940252

This should work:

SELECT a.property_id, a.address, Distance(a.geom, b.geom)
FROM smarter.vow_properties AS a
JOIN (
  SELECT Transform(SetSRID(MakePoint(-74.905868, 39.940252), 4326), 32618) AS geom,
         100 AS dist_m
) AS b ON a.geom && Expand(b.geom, b.dist_m) AND Distance(a.geom, b.geom) <= b.dist_m
ORDER BY Distance;

You could also write the query as you had it.  I sometimes use a
subquery as shown above to avoid repeating the geometry and distance
everywhere I need them.

-- 
Michael Fuhr