How do I perform a proximity search with Postgis?

Question

I downloaded the Geonames database (cities1000) and wrote a small Ruby program to import them into my table (geo_cities). I then added geography column called geog.

Then I converted all of the lat/lon numbers into the geog column using:

update public.geo_cities set geog = st_GeogFromText('SRID=4326;POINT(' || longitude || ' ' || latitude || ')');

Things are looking good. Now what I want to do is find all cities within 100 miles of Prague.

So I can get Prague like:

select * from geo_cities where asciiname = 'Prague' and countrycode = 'CZ';

I'm still learning GIS and Postgres so could someone help me with the simple query?

Answer 1

First, make sure you have an index on your geography column. It will speed up the spatial searches:

CREATE INDEX geo_cities_geog_idx ON geo_cities USING GIST geog;
VACUUM ANALYZE geo_cities(geog);

Then, you can use ST_DWithin (with conversions from miles to metres) on a self-joined query:

SELECT gc.*, ST_Distance(gc.geog, pt.geog)/1609.344 AS distance_miles
FROM geo_cities gc, geo_cities pt
WHERE pt.asciiname = 'Prague' and pt.countrycode = 'CZ'
  AND ST_DWithin(gc.geog, pt.geog, 160934.4)
ORDER BY ST_Distance(gc.geog, pt.geog);

If you are in the market for a good book, check out PostGIS in Action.

Answer 2

For those souls whose DB is with geometry type instead of geography: ST_Transform() or ::geography casting should help.

Example: Suppose we have:

• a location whose 4326 SRID (basic lat-long) value is: POINT(73.8473 18.4869)
• we want to find all the entries within 50 km (50,000 m) of it from a table called table1
• and the table's geometry column is named "geometry". (and it doesn't have a geography column)

Method 1: Using ST_Transform

• the target point is located in India where a good Meters CRS/SRID is 7755

The SQL would be like:

select * from table1 
where ST_DWithin ( 
  ST_Transform(geometry,7755), 
  ST_Transform(ST_GeomFromText('POINT(73.8473 18.4869)',4326) ,7755), 
  50*1000 )
ORDER BY ST_Distance(geometry, ST_GeomFromText('POINT(73.8473 18.4869)',4326));

• So we're applying ST_Transform() to the first 2 arguments of the ST_DWithin() function, to convert them to a CRS/SRID that gives distance in meters.

Method 2: Using ::geography casting

credit to @jbalk for suggesting this.

• Benefit: no need to find out which is the best CRS/SRID for that region.
• target lat-long can directly be cast as geography

select * from table1
where ST_DWithin(
    geometry::geography, 
    ST_GeogFromText('SRID=4326;POINT(73.8473 18.4869)'), 
    50*1000 )
order by ST_Distance(geometry, ST_GeomFromText('POINT(73.8473 18.4869)',4326));

Note on ORDER BY clause:

• No need to apply ST_Transform or ::geography casting in the ORDER BY part,
• as we're just sorting among the candidates that have already qualified for the distance threshold in the WHERE clause
• now the algebraic distances using default lat-long values should compare with each other and give us the same order as what meters conversion would have yielded, at least within a small enough area like a 50km radius.

Refs:

• https://postgis.net/docs/PostGIS_Special_Functions_Index.html
• https://postgis.net/docs/ST_DWithin.html
• https://postgis.net/docs/RT_ST_Transform.html
• https://postgis.net/docs/ST_Distance.html
• https://postgis.net/docs/ST_GeogFromText.html
• https://postgis.net/docs/ST_GeomFromText.html