Community

Blog Events Webinars Tutorials Forum
  1. Blog
  2. Events
  3. Webinars
  4. Tutorials
  5. Forum
Create Account
×
Community Blog PostGIS: How to Calculate Distance between Coordinates in Spherical and Projected Systems

PostGIS: How to Calculate Distance between Coordinates in Spherical and Projected Systems

In this article, the author discusses calculating distances in projected and spherical coordinate systems in respect of geometry and geographic data using PostGIS functions.

By digoal

Background

PostGIS has two popular spatial data types: geometry and geography. What are their differences? Which one should you choose?

In a geographic information system (GIS), two common coordinate system types are used: a global or spherical coordinate system (geographical coordinate system) and a two-dimensional coordinate system (projected coordinate system).

The spherical coordinate system is used mainly for calculation, while the two-dimensional coordinate system is primarily used for demonstration.

When you project the spheroid's surface maps onto a plane, you get the projected coordinates of points on the spheroid's surface. Imagine that the earth is a light orb and wrap it with a cylinder. After projecting Earth's spherical surface maps onto a two-dimensional Cartesian coordinate plane, you get the projected coordinates of the points on Earth's surface on this plane. The larger the projected area, the lower the precision, and vice versa. Of course, the projected sector’s size varies. And different coordinate systems are used in various industries, for example, the coordinate system specific to map drawing.

The most commonly used coordinate systems are the WGS 84 Long Lat (SRID=4326) spherical coordinate system and the Mercator (EPSG:3785) projected coordinate system.

Now, let's look at the geometry and geography types. The geometry type supports both planar and spatial objects. However, the geography type only supports spatial objects.

Geometry supports more functions and has a lower geometric calculation cost.

Geography supports fewer functions and has a higher calculation cost. Then why do we still use it? The reasons are as follows:

4.2.2. When to use Geography Data type over Geometry data type    
    
The geography type allows you to store data in longitude/latitude coordinates,     
but at a cost: there are fewer functions defined on GEOGRAPHY than there are on GEOMETRY;     
those functions that are defined take more CPU time to execute.    
    
The type you choose should be conditioned on the expected working area of the application you are building.     
Will your data span the globe or a large continental area, or is it local to a state, county, or municipality?    
    
If your data is contained in a small area, you might find that choosing an appropriate     
projection and using GEOMETRY is the best solution, in terms of performance and functionality available.    
    
If your data is global or covers a continental region, you may find that GEOGRAPHY     
allows you to build a system without having to worry about projection details.     
You store your data in longitude/latitude, and use the functions that have been defined on GEOGRAPHY.    
    
If you don't understand projections, and you don't want to learn about them,     
and you're prepared to accept the limitations in functionality available in GEOGRAPHY,     
then it might be easier for you to use GEOGRAPHY than GEOMETRY.     
Simply load your data up as longitude/latitude and go from there.      
    
Refer to Section 14.11, “PostGIS Function Support Matrix” for compare between     
what is supported for Geography vs. Geometry.     
For a brief listing and description of Geography functions,     
refer to Section 14.4, “PostGIS Geography Support Functions”

Distance calculation results vary with the projected coordinate system.

If you don't know the projected coordinate system of the geometric points, the distance calculation results are usually inaccurate.

For example, if you transform the following point directly to a 2163 (US National Atlas Equal Area) coordinate system, you get an inaccurate distance.

db1=# SELECT st_distance(ST_Transform(ST_GeomFromText('POINT(120.08 30.96)', 4326), 2163 ), ST_Transform(ST_GeomFromText('POINT(120.08 30.92)', 4326), 2163 ));   
   st_distance          
------------------      
 4030.46766234184      
(1 row)

A 2163 coordinate system looks like this:

postgres=# select * from spatial_ref_sys where srid=2163;    
-[ RECORD 1 ]-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------    
srid      | 2163    
auth_name | EPSG    
auth_srid | 2163    
srtext    | PROJCS["US National Atlas Equal Area",GEOGCS["Unspecified datum based upon the Clarke 1866 Authalic Sphere",DATUM["Not_specified_based_on_Clarke_1866_Authalic_Sphere",SPHEROID["Clarke 1866 Authalic Sphere",6370997,0,AUTHORITY["EPSG","7052"]],AUTHORITY["EPSG","6052"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4052"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",45],PARAMETER["longitude_of_center",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","2163"]]    
proj4text | +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs     
AUTHORITY["EPSG", "9122"]指的是EPSG数据集中UNIT为degree的ID是9122；  
  
AUTHORITY["EPSG", "4326"]指的是地理坐标系WGS 84的ID是4326；  
  
AUTHORITY["EPSG", "9001"]指的是EPSG中UNIT为meter的ID是9001；  
  
AUTHORITY["EPSG", "32650"]指的是该投影坐标系WGS 84 / UTM zone 50N的ID是32650

This presents a false impression:

用st_distance函数计算出来的经纬度之间的距离，跟用java程序算出来的距离相差很大。      
      
这个例子，st_distance算出来的距离是4030，我们程序算出来的是4445，换另外两个相距很远的点，这个差距值会更大。

Solutions

The answer is simple: calculate the spherical distance instead of the planar distance.

1) Using ST_DistanceSpheroid to Calculate the Spherical Distance Between Two Geometric Points in a Spherical Coordinate System

You can use ST_DistanceSpheroid to calculate the distance between two geometric points based on their spherical coordinates. The spheroid attributes are automatically configured.

SPHEROID["Krasovsky_1940",6378245.000000,298.299997264589]

postgres=# SELECT st_distancespheroid(ST_GeomFromText('POINT(120.08 30.96)', 4326),ST_GeomFromText('POINT(120.08 30.92)', 4326), 'SPHEROID["WGS84",6378137,298.257223563]');    
 st_distancespheroid     
---------------------    
    4434.73734584354    
(1 row)

2) Using ST_Distance to Calculate the Distance Between Two Geometric Points on the Projected Plane in a Projected Coordinate System

To ensure that the results are accurate, you must use accurate projected coordinates in a small projected coordinate system that is large enough to cover both the points to be calculated.

ST_Transform(ST_GeomFromText('POINT(120.08 30.96)', 4326), 2163 )

If the geometric value’s spatial reference system identifier (SRID) does not match the SRID of the (high-accuracy) desired coordinate system, you can use the ST_Transform function to convert it into desired coordinates before calculating the distance.

postgres=# SELECT st_distance(ST_Transform(ST_GeomFromText('POINT(120.08 30.96)', 4326), 2390 ), ST_Transform(ST_GeomFromText('POINT(120.08 30.92)', 4326), 2390 ));   
   st_distance      
------------------  
 4547.55477647394  
(1 row)

3) Using ST_Distance to Calculate the Spherical Distance of Two Geographic Points in a Spherical Coordinate System

float ST_Distance(geography gg1, geography gg2, boolean use_spheroid);   -- 适用椭球体（WGS84）    
    
    
use_spheroid设置为ture表示使用:      
  -- WGS84 椭球体参数定义    
  vspheroid := 'SPHEROID["WGS84",6378137,298.257223563]' ;

In the following snippet, replace "XXXX" with the SRID of your spherical coordinate system. For more information about available coordinate systems, see the spatial_ref_sys table.

postgres=# SELECT st_distance(ST_GeogFromText('SRID=xxxx;POINT(120.08 30.96)'), ST_GeogFromText('SRID=xxxx;POINT(120.08 30.92)'), true);    
  st_distance       
----------------    
 xxxxxxxxxxxxxxxxxxxx    
(1 row)    
  
例如  
  
postgres=# SELECT st_distance(ST_GeogFromText('SRID=4610;POINT(120.08 30.96)'), ST_GeogFromText('SRID=4610;POINT(120.08 30.92)'), true);    
  st_distance     
----------------  
 4434.739418211  
(1 row)

If some deviation is allowed, you can use the 4326 (WGS 84) spherical coordinate system.

postgres=# SELECT st_distance(ST_GeogFromText('SRID=4326;POINT(120.08 30.96)'), ST_GeogFromText('SRID=4326;POINT(120.08 30.92)'), true);    
  st_distance       
----------------    
 4434.737345844    
(1 row)

The geography type supports only spherical coordinate systems, so an error will be returned if a projected coordinate system is used.

postgres=# SELECT st_distance(ST_GeogFromText('SRID=2369;POINT(120.08 30.96)'), ST_GeogFromText('SRID=2369;POINT(120.08 30.92)'), true);    
错误:  22023: Only lon/lat coordinate systems are supported in geography.  
LOCATION:  srid_is_latlong, lwgeom_transform.c:774

4) Using ST_DistanceSpheroid to Calculate the Spherical Distance Between Two Geographic Points in a Spherical Coordinate System by Specifying the Spheroid Parameters

This approach provides the most accurate result, but you must know the landscape’s geographical properties. Specifically, you must specify the spheroid parameters.

db1=# SELECT st_distancespheroid(ST_GeomFromText('POINT(120.08 30.96)', 4326),ST_GeomFromText('POINT(120.08 30.92)', 4326), 'SPHEROID["WGS84",6378137,298.257223563]');        
 st_distancesphere       
-------------------      
    4447.803189385      
(1 row)

Summary

When calculating the distance between two points, consider that they're on a spheroid to get an accurate result.

When choosing the data type (geometry or geography), we recommend using geometry, which supports both the spherical and the projected coordinate systems. You can choose the most suitable coordinate system based on your knowledge of the points’ geographical locations.

-- 适用平面直角坐标，适用geometry类型，计算直线距离    
float ST_Distance(geometry g1, geometry g2);             
    
-- 适用椭球体坐标，适用geography类型，计算球面距离  
float ST_Distance(geography gg1, geography gg2);         
    
-- 适用椭球体坐标（WGS84），适用geography类型，计算球面距离    
float ST_Distance(geography gg1, geography gg2, boolean use_spheroid);     
    
   use_spheroid设置为ture表示使用:      
     vspheroid := 'SPHEROID["WGS84",6378137,298.257223563]' ; -- WGS84 椭球体参数定义     
    
-- 适用椭球体坐标以及投影坐标，适用geometry类型，求球面距离（需要输入spheroid）  
float ST_DistanceSpheroid(geometry geomlonlatA, geometry geomlonlatB, spheroid measurement_spheroid);

References

1) ST_DistanceSphere: Calculation of Spherical Distances

2) ST_Distance: Calculation of Spherical or Planar Distances (Depending on the Input Type)

3) Transformation of Coordinate Systems

4) http://blog.163.com/jey_df/blog/static/18255016120149145755781/

5) https://desktop.arcgis.com/zh-cn/arcmap/10.3/guide-books/map-projections/mercator.htm

PostgreSQL Database PostGIS Spherical Coordinates Planar Coordinates Spherical Distance Planar Distance
0 0 0
Share on

Read previous post:

PostgreSQL: Importing OpenStreetMap and Road Data

Read next post:

PostgreSQL 12: B-tree Index Improvements – Duplicate Key and Sort by CTID

digoal

281 posts | 24 followers

You may also like

Comments

digoal

281 posts | 24 followers

Related Products

More Posts by digoal

See All