By Digoal
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
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)
PostGIS Long Lat Geometry Distance Search Tuning Using the GiST KNN Function
digoal - February 16, 2021
digoal - December 14, 2020
digoal - December 18, 2020
digoal - December 23, 2020
digoal - January 25, 2021
digoal - May 16, 2019
Alibaba Cloud PolarDB for PostgreSQL is an in-house relational database service 100% compatible with PostgreSQL and highly compatible with the Oracle syntax.
Learn MoreAn online MPP warehousing service based on the Greenplum Database open source program
Learn MoreAn on-demand database hosting service for PostgreSQL with automated monitoring, backup and disaster recovery capabilities
Learn MoreLeverage cloud-native database solutions dedicated for FinTech.
Learn MoreMore Posts by digoal