Calculating accurate distances using various ST_Distance options

Question

I want to get the nearest neighbour information for a point, and I'm getting different results based on different geometry/geography types. Assuming I only care about the most accurate result, which version should I use?

1. Using a sphere: SELECT ('SRID=4326;POINT(-95.1212742 29.5574348)'::geography <-> 'SRID=4326;POINT(-91.0717865 30.3703801)'::geography)/1000 --RETURNS 400.40
2. Using a planar distance (not sure what that means for ESPG:4326): SELECT ST_Distance('SRID=4326;POINT(-95.1212742 29.5574348)'::geometry, 'SRID=4326;POINT(-91.0717865 30.3703801)'::geometry) --RETURNS 4.130 <- this is identical to using <-> operator on geometry type.
3. Using a projected coordinate system: SELECT ST_Distance('SRID=4326;POINT(-95.1212742 29.5574348)'::geography, 'SRID=4326;POINT(-91.0717865 30.3703801)'::geography)/1000 --RETURNS 401.09 <- geodesic distance based on SRID

We can ignore #1 because it uses the fastest spherical calculation. However, for the last two, I get diverging results given different neighbour points, and I'm not sure which one is the more accurate one:

SELECT ST_Distance('SRID=4326;POINT(-95.1212742 29.5574348)'::geometry, 'SRID=4326;POINT(-91.1505645 30.6506557)'::geometry) --RETURNS 4.118
SELECT ST_Distance('SRID=4326;POINT(-95.1212742 29.5574348)'::geography, 'SRID=4326;POINT(-91.1505645 30.6506557)'::geography)/1000 --RETURNS 401.42

SELECT ST_Distance('SRID=4326;POINT(-95.1212742 29.5574348)'::geometry, 'SRID=4326;POINT(-91.0717865 30.3703801)'::geometry) --RETURNS 4.130
SELECT ST_Distance('SRID=4326;POINT(-95.1212742 29.5574348)'::geography, 'SRID=4326;POINT(-91.0717865 30.3703801)'::geography)/1000 --RETURNS 401.09

#2 Returns 2D Cartesian (planar) distance between two geometries, in projected units. Considering EPSG:4326 is not a planar CRS, this will return result in degrees, which would be less accurate?

#3 returns geodesic distance based on SRID, which should be most accurate?

Answer 1

Accuracy in descending order:

1. Spheroidal distance - calculated on the reference ellipsoid of the given SRID:

   SELECT
     ST_Distance(
       'SRID=4326;POINT(-95.1212742 29.5574348)'::GEOGRAPHY, 
       'SRID=4326;POINT(-91.0717865 30.3703801)'::GEOGRAPHY,
       use_spheroid = TRUE  -- DEFAULT
     )
   ;
   

2. 1. Planar distance (suitably projected) - calculated on a Cartesian plane using a specific projection on a large scale:

      SELECT
        ST_Distance(
          ST_Transform(
            'SRID=4326;POINT(-95.1212742 29.5574348)'::GEOMETRY,
            <SRID>
          ),
          ST_Transform(
            'SRID=4326;POINT(-91.0717865 30.3703801)'::GEOMETRY,
            <SRID>
          )
        )
      ;
      
      SELECT
        ST_Transform('SRID=4326;POINT(-95.1212742 29.5574348)'::GEOMETRY, <SRID>) <-> ST_Transform('SRID=4326;POINT(-91.0717865 30.3703801)'::GEOMETRY, <SRID>)
      ;
      

   2. Spherical distance - calculated on a sphere:

      SELECT
        ST_Distance(
          'SRID=4326;POINT(-95.1212742 29.5574348)'::GEOGRAPHY, 
          'SRID=4326;POINT(-91.0717865 30.3703801)'::GEOGRAPHY,
          use_spheroid = FALSE
        )
      ;
      
      SELECT
        'SRID=4326;POINT(-95.1212742 29.5574348)'::GEOGRAPHY <-> 'SRID=4326;POINT(-91.0717865 30.3703801)'::GEOGRAPHY
      ;
      

3. Planar distance (unprojected, or unsuitably projected) - calculated on a Cartesian plane:

   SELECT
     ST_Distance(
       'SRID=4326;POINT(-95.1212742 29.5574348)'::GEOMETRY, 
       'SRID=4326;POINT(-91.0717865 30.3703801)'::GEOMETRY
     )
   ;
   
   SELECT
     'SRID=4326;POINT(-95.1212742 29.5574348)'::GEOMETRY <-> 'SRID=4326;POINT(-91.0717865 30.3703801)'::GEOMETRY
   ;
   

   Due to the cosinuidal degradion of the relation between longitudinal and latitudinal surface distance per increasing latitude, geographic distances (in degree) are meaningless without its angular context.

Answer 2

The most accurate would be #3:

SELECT 
  ST_Distance(
    'SRID=4326;POINT(-95.1212742 29.5574348)'::geography, 
    'SRID=4326;POINT(-91.0717865 30.3703801)'::geography
  )

#1 is the operator to do knn search, and yes it uses sphere computation for geography

#2 never use geometry for 4326, you can do that only if you project your data first. The result you get if you use it as geometry is in degree, which means nothing because depending where you are on the globe lat and lon degrees don't have the same size, so mixing them as if they were equivalent is an error. If you project your data first, your knn search would be faster too, and the precision would depends of the accuracy of the projection for your working zone.