Boundary Issues in PostGIS Coordinate Transformation (SRID) - ST_Transform

By Digoal

Background

When using PostgreSQL ST_Transform to convert coordinates, a user encountered a boundary issue, causing an inaccurate distance calculation.

Sheyu from Cainiao Network provided a perfect solution for this issue.

It is not about whether to use the 26986 coordinate system. First, understand how the 26986 coordinate system is defined. Why does it function well when calculating a negative longitude (Western Hemisphere)? This is because 26986 is Massachusetts' local coordinate system (regional coordinate system), and its two main definition parameter origins are the longitude and latitude (-71.5, 41). Thus, the central meridian is 71.5°W. It is far away from China (China is centered on 120°E), and therefore the projected coordinate system is inapplicable. Simply put, the farther away from the central meridian, the greater the projection error will be.

To calculate the distance between latitude and longitude on a spherical surface, use functions provided by PostGIS without projecting onto a plane first, unless you want to perform other plane operations such as area calculation. In this case, how to select the projection SRID?

Note that a geographic coordinate system and a projected coordinate system are different. The former is based on a spheroid, and the latter is based on a plane.

Geographic coordinate system:

4214 GCS_Beijing_1954  
4326 GCS_WGS_1984   
4490 GCS_China_Geodetic_Coordinate_System_2000  
4555 GCS_New_Beijing  
4610 GCS_Xian_1980

Projected coordinate system:

2327 Xian_1980_GK_Zone_13 
2328 Xian_1980_GK_Zone_14 
2329 Xian_1980_GK_Zone_15 
2330 Xian_1980_GK_Zone_16 
2331 Xian_1980_GK_Zone_17 
2332 Xian_1980_GK_Zone_18 
2333 Xian_1980_GK_Zone_19 
2334 Xian_1980_GK_Zone_20 
2335 Xian_1980_GK_Zone_21
2336 Xian_1980_GK_Zone_22 
2337 Xian_1980_GK_Zone_23 
2338 Xian_1980_GK_CM_75E 
2339 Xian_1980_GK_CM_81E 
2340 Xian_1980_GK_CM_87E 
2341 Xian_1980_GK_CM_93E 
2342 Xian_1980_GK_CM_99E 
2343 Xian_1980_GK_CM_105E 
2344 Xian_1980_GK_CM_111E 
2345 Xian_1980_GK_CM_117E 
2346 Xian_1980_GK_CM_123E 
2347 Xian_1980_GK_CM_129E 
2348 Xian_1980_GK_CM_135E  
2349 Xian_1980_3_Degree_GK_Zone_25 
2350 Xian_1980_3_Degree_GK_Zone_26 
2351 Xian_1980_3_Degree_GK_Zone_27 
2352 Xian_1980_3_Degree_GK_Zone_28 
2353 Xian_1980_3_Degree_GK_Zone_29 
2354 Xian_1980_3_Degree_GK_Zone_30 
2355 Xian_1980_3_Degree_GK_Zone_31 
2356 Xian_1980_3_Degree_GK_Zone_32 
2357 Xian_1980_3_Degree_GK_Zone_33 
2358 Xian_1980_3_Degree_GK_Zone_34 
2359 Xian_1980_3_Degree_GK_Zone_35 
2360 Xian_1980_3_Degree_GK_Zone_36 
2361 Xian_1980_3_Degree_GK_Zone_37 
2362 Xian_1980_3_Degree_GK_Zone_38 
2363 Xian_1980_3_Degree_GK_Zone_39 
2364 Xian_1980_3_Degree_GK_Zone_40 
2365 Xian_1980_3_Degree_GK_Zone_41 
2366 Xian_1980_3_Degree_GK_Zone_42 
2367 Xian_1980_3_Degree_GK_Zone_43 
2368 Xian_1980_3_Degree_GK_Zone_44 
2369 Xian_1980_3_Degree_GK_Zone_45

Example

The following two 4326 coordinates are slightly different. But when converted into 26986 coordinates, their difference becomes significant.

postgres=# select ST_AsEWKT(ST_Transform(ST_GeomFromText('POINT(108.50000000001 22.8)', 4326), 26986));  
                      st_asewkt                        
-----------------------------------------------------  
 SRID=26986;POINT(8123333.59043839 12671815.6459695)  
(1 row)  
  
postgres=# select ST_AsEWKT(ST_Transform(ST_GeomFromText('POINT(108.5000000001 22.8)', 4326), 26986));  
                      st_asewkt                         
------------------------------------------------------  
 SRID=26986;POINT(-7723333.59044452 12671815.6459593)  
(1 row)

The distance calculated based on the converted coordinates is inaccurate because the central meridian of the SRID is far away from the calculation point. The projection error increases when the distance from the central meridian increases.

postgres=# select * from spatial_ref_sys where srid=26986;
-[ RECORD 1 ]-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
srid      | 26986
auth_name | EPSG
auth_srid | 26986
srtext    | PROJCS["NAD83 / Massachusetts Mainland",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",42.68333333333333],PARAMETER["standard_parallel_2",41.71666666666667],PARAMETER["latitude_of_origin",41],PARAMETER["central_meridian",-71.5],PARAMETER["false_easting",200000],PARAMETER["false_northing",750000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","26986"]]
proj4text | +proj=lcc +lat_1=42.68333333333333 +lat_2=41.71666666666667 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +datum=NAD83 +units=m +no_defs 


postgres=# select * from spatial_ref_sys where srid=4326;
-[ RECORD 1 ]---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
srid      | 4326
auth_name | EPSG
auth_srid | 4326
srtext    | GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
proj4text | +proj=longlat +datum=WGS84 +no_defs

PROJCS indicates projected coordinates, and GEOGCS indicates spherical coordinates.

To solve this problem, do not convert the coordinates into 26986 coordinates (planar coordinates). Use the corresponding function of PostGIS to calculate the spherical coordinates. Even if you need to calculate a planar distance, use a coordinate system commonly used in China (whose central meridian is in China, as described earlier).

http://postgis.net/docs/manual-2.0/ST_Distance.html

try this:  
  
postgres=# select ST_Distance(ST_GeographyFromText('SRID=4326;POINT(108.51 22.8)'), ST_GeographyFromText('SRID=4326;POINT(108.499999999999999 22.79)'));  
-[ RECORD 1 ]--------------------  
st_distance | 1510.16913796499989  
  
-- Geography example -- same but note units in meters - use sphere for slightly faster less accurate  
  
-- Geometry example - units in meters (SRID: 26986 Massachusetts state plane meters) (most accurate for Massachusetts)

Community

Boundary Issues in PostGIS Coordinate Transformation (SRID) - ST_Transform

Background

Example

Read previous post:

Read next post:

digoal

You may also like

Comments

digoal

Related Products

PolarDB for PostgreSQL

Database Overview

AnalyticDB for PostgreSQL

ApsaraDB RDS for PostgreSQL