How to calculate distance between two 3D points in kilometers using postgis

Question

I get latitude, longitude and altitude using android.location.Location, after that I send this data on the server, and save it using:

update users set location = ST_SetSRID(ST_MakePoint(-40, 40, 40), 4326) where id = 1;

users.location column was created using:

SELECT AddGeometryColumn ('users', 'location', 4326, 'POINT', 3);

I already saw an answers about 2D points, or answers in which you need to translate one SRID into another, and this is confused me. I want to know correct way from the Android to the postgis.

So, how can I calculate distance between user with id 1 and user with id 2 in kilometers using postgis?

UPD:

SELECT AddGeometryColumn ('users', 'location', 4326, 'POINT', 2);

update users set location = ST_SetSRID(ST_MakePoint(-118.4079, 33.9434), 4326) where id = 1;

update users set location = ST_SetSRID(ST_MakePoint(2.5559, 49.0083), 4326) where id = 2;

SELECT ST_Distance(a.location, b.location)/1000 as Dist_deg from users a, users b where a.id=1 AND b.id=2;

Output:

 0.121898285970107

This is wrong(since these are degrees), so I tried this:

SELECT ST_Distance(ST_Transform(a.location, 2163), ST_Transform(b.location, 2163))/1000 as Dist_deg from users a, users b where a.id=1 AND b.id=2;

Output:

 8731.12032614592

So, this is wrong too.

Answer 1

For the distance between the 3D points, you should create the 3d geometry for these points, then transform the points on an SRID (depends on the location of the points). After that, you can use (ST_3DDISTANCE), which is the only way.

You can find this ANSWER for a similar problem as yours.

If you want to drop the Z-Index, you can calculate the distance between the two 2D points, but be careful that the calculated distance (in degrees which is nonsense number) because the units for spatial reference 4326 are degrees.

 SELECT ST_Distance(
      ST_GeometryFromText('POINT(-118.4079 33.9434)', 4326), -- Los Angeles (LAX)
      ST_GeometryFromText('POINT(2.5559 49.0083)', 4326)     -- Paris (CDG)
      ) as Dist_deg;

In order to calculate the distance (in meter or Km), we must treat geographic coordinates not as approximate Cartesian coordinates but rather as true spherical coordinates. We must measure the distances between points as true paths over a sphere – a portion of a great circle [1].

So in this case, we will use (ST_Geography) as follows

 SELECT ST_Distance(
  ST_GeographyFromText('POINT(-118.4079 33.9434)'), -- Los Angeles (LAX)
  ST_GeographyFromText('POINT(2.5559 49.0083)')     -- Paris (CDG)
  )/1000 as Dist_KM;

And the result is about 9124.665 Km

For your case e.g (Location between 2D point1 & 2D point2):

SELECT ST_Distance(a.location, b.location)/1000 as Dist_deg from users a, users b where a.id=1 AND b.id=2;

For more information [please find this helpful tutorial]

Finally, let's have a small comparison as an example between the 3D & 2D distance[2].

SELECT ST_3DDistance(
            ST_Transform(ST_GeomFromEWKT('SRID=4326;POINT(-72.1235 42.3521 4)'),2163),
            ST_Transform(ST_GeomFromEWKT('SRID=4326;LINESTRING(-72.1260 42.45 15, -72.123 42.1546 20)'),2163)
        ) As dist_3d,
        ST_Distance(
            ST_Transform(ST_GeomFromText('POINT(-72.1235 42.3521)',4326),2163),
            ST_Transform(ST_GeomFromText('LINESTRING(-72.1260 42.45, -72.123 42.1546)', 4326),2163)
        ) As dist_2d;

And the results are (3D Distance=127.295059324629, 2D Distance=126.66425605671)

Answer 2

Much depends on what scale you are working.

• For distances on global scale, the ad hoc measurements for true 3D space in PostGIS are limited by projections and their distortions; ST_3DDistance needs all coordinates projected in meter to make sense, and accuracy of worldwide projections is very modest on global scale (check e.g. this answer for possible options).
  Considering the overall bad accuracy for mobile based GPS, especially for elevation, on that scale a 'simple' spheroid based 2D distance would be more accurate than anything involving projections I guess.
• On larger scales however, a simple transformation into the appropriate UTM projection could be all you need (the API docs on android.location.getAltitude states elevation is in meter above WGS84 ellipsoid, thus having a WGS84 based projection makes sense). There are other, local projections with possibly lesser distortion, but if your users are distributed around the world (not to be confused with a measurement around the globe; this here means two user to be measured are in the same area somewhere on Earth), there is e.g. this function to find the proper UTM zone and SRID for given WGS84 LatLons, making that a lot easier to manage (possibly at the cost of slightly less accuracy). GPS elevation is still no better that an educated guess, but here ST_3DDistance could work well.

Assuming the latter, to make your data usable you could

• transform points on the fly; possibly slowest in execution for many measurements
• update the geometry column or create a second one; possible drawback that you have to run the update on every insert
• use triggers; do the above on insert

following this general workflow:

• creating 2D Point from LatLon
• (finding UTM Zone and) transforming them to proper (UTM) projection
• create POINT Z from POINTs with elevation assigned
• 3DDistance - energize!

For 2D distances (with geography type), the same options arise, but measuring is straight forward.

For a 3D distance measurement between Los Angeles and San Diego on UTM projection (WGS 84 / UTM zone 11S - EPSG:32711), running (thanks to @User2009)

WITH
    pts AS (
        SELECT ST_GeographyFromText('POINT(-118.4079 33.9434)') AS a_geom_wgs, -- Los Angeles (LAX)
               ST_GeographyFromText('POINT(-117.1625 32.715)') AS b_geom_wgs,    -- San Diego
               ST_Transform(ST_GeometryFromText('POINT(-118.4079 33.9434)', 4326), 32611)  AS a_geom_utm, -- Los Angeles (LAX)
               ST_Transform(ST_GeometryFromText('POINT(-117.1625 32.715)', 4326), 32611) AS b_geom_utm    -- San Diego
    )
SELECT ST_Distance(a_geom_wgs::geography, b_geom_wgs::geography) / 1000 AS dist_spheroid_2d,
       ST_3DDistance(ST_MakePoint(ST_X(a_geom_utm), ST_Y(a_geom_utm), 0), ST_MakePoint(ST_X(b_geom_utm), ST_Y(b_geom_utm), 0)) / 1000 AS dist_utm_3d
FROM pts

with elevation set to

• 0m for both points returns dist_utm_3d = 178.843582335076 (and dist_spheroid_2d = 178.90107074162 as a comparison),
• 3000m for Los Angeles returns dist_utm_3d = 178.868742217424

...which is a little disappointing, but does it´s job (using the same place but different altitudes returns the altitude as the correct distance).

Choose your weapon.