Search code examples
postgresqlpostgis

PostGIS: how to convert a point column into a type useable by ST_DWithin?


I have the following query to try and get a total of places within the range of a provided point, comparing against my table of places. All places have a geo column, defined as geo point default null.

select
  count(*)
from
  places
where
  ST_DWithin(
    ST_GeomFromText(concat('POINT', geo::text), 4326),
    ST_GeomFromText('POINT(37.64903402157866, -83.84765625000001)', 4326),
    66 * 1609.34
  );

However, that returns:

parse error - invalid geometry

I tried casting to ::geometry, using other functions, but nothing really works. Here's how my data looks like:

select geo, geo::text, geo::geometry from places where geo is not null limit 10

returns:

{"x":-84.3871,"y":33.7485}
(-84.3871,33.7485)
01010000000612143FC61855C02B8716D9CEDF4040

Am I doing something wrong? Should I alter the table and store another column type? How do I not lose the existing data?


Solution

    1. Your lon/lat order is flipped between what you're looking for and what you have in the table. You'll need to make sure which order is correct and flip back the target point or the data in the table. PostGIS operates in x+y, lon+lat, while Google Maps is an example where you can find y+x, lat+lon. Always check what you're dealing with.

    2. 66 * 1609.34 seems like conversion from miles to meters, but since you're using type geometry in 4326, STDWithin() will assume a degree to be your unit of measurement:

      For geometry: The distance is specified in units defined by the spatial reference system of the geometries.

      Which means your target range exceeds the circumference of the globe, around 295 times.

    3. Postgres' built-in point is not the same type as geometry(Point) from PostGIS extension. You'll need to cast or convert between them.

    4. STDWithin() can measure in meters, but only if you give it geography types (or switch the CRS/SRID to meter-based), so that's another cast. It can also increase the precision at the cost of some performance, if you tell it to assume a spheroid, as it by default assumes the globe is a regular sphere.

    Demo at db<>fiddle:

    create table places(id int generated by default as identity primary key, 
                        geo point,
                        comment text);
    insert into places values 
      (1,point(-84.3871,33.7485),'your existing point'),
      (2,point(33.7485,-84.3871),'flipped lat/lon of existing point'),
      (3,point(37.64903402157866, -83.84765625000001),'target coords'),
      (4,point(-83.84765625000001, 37.64903402157866),'flipped target coords');
    
    select
      count(*),array_agg(id::text||': '||comment)
    from
      places
    where
      ST_DWithin(
        geo::geometry::geography(Point,4326),
        ST_GeogFromText('SRID=4326;POINT(37.64903402157866 -83.84765625000001)'),
        66 * 1609.34, --66 miles
        true --measure over spheroid
      );
    
    count array_agg
    2 {"2: flipped lat/lon of existing point",
    "3: target coords"}