Postgres returning in-correct results with distance join queries

Question

I created a table containing locations along with their spatial coordinates expressed in (latitude, longitude (in degrees)) in postgres. Commands which I used for the same are:

I inserted spatial coordinates as (longitude, latitude) in postgres

I created a table containing locations along with their spatial coordinates expressed in (latitude, longitude (in degrees)) in postgres. Commands which I used for the same are:

CREATE TABLE spatialTest(
  name character varying(100),
  the_geo geography
);
\copy spatialTest(name,the_geo) FROM 'testSpatial.csv' DELIMITERS E'\t' CSV;

testSpatial.csv contains the following values:

A   SRID=4326;POLYGON((0.178773 -127.681841,0.178711 -127.681962,0.179125 -127.682083,0.179176 -127.682006,0.179153 -127.681986,0.179143 -127.681962,0.179147 -127.681935,0.179166 -127.681913,0.179195 -127.681897,0.179244 -127.681886,0.179284 -127.681887,0.179336 -127.681904,0.179464 -127.681757,0.179489 -127.681736,0.179429 -127.681680,0.179370 -127.681516,0.179221 -127.681331,0.179184 -127.681185,0.179051 -127.681264,0.178822 -127.681499,0.178761 -127.681698,0.178796 -127.681703,0.178839 -127.681721,0.178857 -127.681736,0.178861 -127.681740,0.178871 -127.681756,0.178873 -127.681782,0.178859 -127.681809,0.178843 -127.681825,0.178812 -127.681839,0.178773 -127.681841))
B   SRID=4326;POINT(0.628912 -127.700922)

Now I want to find all spatial locations which are within a distance of 50 km of each other. For doing so I used the following command:

select s1.name, s2.name  from spatialTest s1, 
     spatialTest s2  where ST_DWithin(s1.the_geo, s2.the_geo, 50000);

However, to my surprise I found that although A and B are separated from each other by a distance greater than 50 km (50.0995 km to be precise. Found using Chris Veness's geodesy formulae (Calculate distance between a point and a line segment in latitude and longitude)), yet they are returned by postgres as results. Can someone please help me figure out as to where am I going wrong.

I am using PostgreSQL 9.6devel and Postgis version which I am using is: POSTGIS="2.2.1 r14555"

You got an answer for that on the Postgres "Bugs" mailing list: "I think you need to ask about that on the Postgis mailing lists, not here" - did you do that? And why aren't you using a proper Postgres release? — a_horse_with_no_name
Your code and sample data return a error by me: ERROR: parse error - invalid geometry HINT: "SRID=4326;POINT(0.628912,-" <-- parse error at position 26 within geometry. Please add correct data and code examples — Tom-db
"spatial coordinates expressed in (latitude, longitude (in degrees)) in postgres" I think you need to pass arguments in longitude, latitude order. That is, first arg is longitude, not latitude. Could this be all or part of your problem? — Mike Sherrill 'Cat Recall'
you should post it on the PostGIS mailing list, not the Postgres mailing list — a_horse_with_no_name
Your coordinates are invalid lon/lat coordinates. Maybe are they in the reversed order? COPY returns this warning: NOTICE: Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY — Tom-db

Evan Carroll Evan Carroll · Accepted Answer · 2017-01-06T08:14:33

As others have said, you're getting this on load..

NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY
CONTEXT:  COPY spatialtest, line 2, column the_geo: "SRID=4326;POINT(0.628912 -127.700922)"

Likely, you're trying to load (lat,long), and PostgreSQL expects (long,lat). However we can work around the ordering of your input column. All we have to do is,

Create a temp table asgeometry`.

CREATE TEMP TABLE foo(name text, the_geo geometry);

Chop off the SRID bit (EKWT) losing the SRID on import.
```
sed -i -e's/SRID=4326;//' ./foo.csv
```

Load the input.

\copy foo(name,the_geo) FROM 'foo.csv' DELIMITERS E'\t' CSV;

Reorder the points.

Copy into main table.

INSERT INTO spatialTest(name, the_geo)
SELECT
  name,
  ST_Union(
    ST_SetSRID(
      ST_MakePoint(st_y(geom), st_x(geom))
    ,4326)
  )::geography
FROM foo
CROSS JOIN ST_DumpPoints(foo.the_geo) AS t
GROUP BY name
ORDER BY name;

As a side note: don't use varchar for name. Just leave it as TEXT. Unless you have a good reason to impose a character limit.

Then run your distance..

SELECT ST_Distance(
  (SELECT the_geo FROM spatialtest WHERE name='A'), (SELECT the_geo FROM spatialtest WHERE name='B')
);
  st_distance   
----------------
 49740.52214618
(1 row)

Postgres returning in-correct results with distance join queries

1 Answers