<#> doesn't give a useful distance [duplicate]

Question

I'm new to postgresql/postgis and trying to measure a distance. So st_distance() may be the most accurate, it's the slowest too (as I've read). I read that one can use geometric (instead of geographic) measuring with less accuracy, but the numbers I get don't make sense. This is what I've got:

(all points are geography(Point,4326))

select
  st_astext(point),
  st_astext(point2),
  round(st_distance(point, point2)),
  round(point::geometry <#> point2::geometry) -- cast to geometry and use distance operator
from testgeo2 limit 50

"POINT(53.887 54.043)","POINT(70.764 80.922)","3053406","32" -- ok
"POINT(93.262 59.546)","POINT(84.331 84.334)","2774883","26" -- ok
"POINT(99.624 67.098)","POINT(95.718 85.252)","2027908","19" -- ok
"POINT(54.406 84.531)","POINT(80.27  83.385)", "326308","26" -- hm?!

So I don't get how this can be right. Any suggestions?

One more example

-- new york, tokyo
insert into testgeo2(point, point2) values (
   st_geographyfromtext('point(-74.007873 40.717602)'), 
   st_geographyfromtext('point(139.754333 35.675853)'))

-- get distance
select
  st_astext(point),
  st_astext(point2),
  -- km, most accuracy
  round(st_distance(point, point2)) / 1000,
  -- km, with less accuracy
  round(st_distance(point, point2, false)) / 1000,
  -- km, Lambert
  round(st_distance(
    st_transform(point::geometry, 3587),
    st_transform(point2::geometry, 3587))) / 1000,
  -- km, from degree
  round(point::geometry <#> point2::geometry) * 111
from testgeo2

-- results
POINT(-74.007873 40.717602) POINT(139.754333 35.675853), 
   10871, 10847, 11224, 23754

is geometry just not working on a larger scale?

Answer 1

Considering your first example, you do not appear to be comparing like with like. See what happens when you make it do more:

select
st_astext (point),
st_astext (point2),
round (st_distance (point, point2)),
round (st_distance (point::geometry, point2::geometry)),
round (point::geometry <#> point2::geometry)
from testgeo2 limit 5

Note that ST_Distance (geog, geog) returns "spheroidal minimum distance between two geographies in meters". However, ST_Distance (geom, geom) returns plane minimum distance between two geometries "in projected units", which in your case, will be degrees. And point <#> point is "almost the same as distance" (pnt, pnt).

If you "just" cast a geography to a geometry all you're doing is treating geographic coords as though they were planar and getting the Euclidean, planar distance.

Now for your 2nd example. The two variations of ST_Distance (geog, geog), respectively, return 1) the (default) ellipsoidal, and as you correctly state, most accurate, distance, and 2) the faster, spherical distance. And they agree to within 24km.

For the next distance, you transform the points onto Michigan's "state plane" system via Lambert's conical conformal projection, with two standard parallels at 44.2 and 45.7 degrees (epsg:3587). Now, only along those two parallels will the scale be true and yet the great circle route from new york to tokyo (assuming that's where the points are) veers from 40.7 to 35.7 degrees of latitude. So, you can expect the distance to be in error, unless you determine and apply the appropriate scale factor (see calculating-distance-scale-factor-by-latitude-for-mercator for a related example).

Finally, the 4th distance yields degrees, which you convert to km via 111km/degree. However, that "scale factor" only applies to all north-south distances (and along the equator). You'd really need to calculate, say, three such scale factors -- at beginning, middle and end of the route -- and use a weighted average.

Your suggestion that calculating distances between "geometries" just does not work (without distortion) over larger distances is, generally, and often spectacularly true! That, and the fact that the <#> operator doesn't give geographic distances, is not the fault of PostGIS; it's a matter of understanding the principles and limitations of map projections and the practicalities of planar versus spheroidal COGO within spatial databases.