11

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?

2 Answers 2

13

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.

1

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:

2
  • You could also cast the geometry to geography st_distance(geom1::geography, geom2::geography)
    – jbalk
    Commented Jul 13, 2022 at 18:15
  • @jbalk nice! this looks more elegant, plus we don't have to worry about finding the right meters SRID. Changing my answer. Thanks!
    – Nikhil VJ
    Commented Jul 14, 2022 at 4:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.