作者
digoal
日期
2017-06-22
标签
PostgreSQL , PostGIS , ST_Transform , SRID , 26986
背景
某个用户在使用PostgreSQL ST_Transform转换坐标时，遇到一个边界问题（暂时不清楚是不是BUG，因为对SRID还不太了解），导致距离计算不准确。
菜鸟的舍余同学给出了专业的解答：
不是用不用26986的问题，可以先了解下26986坐标系是怎么定义的哈，为什么它计算负经度（西半球）时没问题呢？这是因为26986为美国马萨诸塞州地方坐标系（区域坐标系），其两个主要定义参数原点为经纬度（-71.5，41），也就是中央子午线为-71.5度，离中国太远了（中国是以东经120度为中心的），不能随便用投影坐标系呀（简单解释就是，离中央子午线越远的地方投影误差越大）
如果是计算球面经纬度距离，postgis有现成函数，无需先投影成平面哈，除非进行其他平面操作，如计算面积，那怎么选投影srid呢？参考https://wenku.baidu.com/view/5975ff2ebe23482fb5da4c24.html， 这里需要注意下地理坐标系和投影坐标系的区别（一个球面一个平面的）。
更多背景知识请参考：
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
例子
下面两个4326坐标系的坐标，只相差了一点点，但是转换为26986坐标系时，出现了翻转。
```
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)
```
使用转换后的坐标，计算距离，导致数据不准确，原因就是这个SRID的中央子午线离计算点太远，离中央子午线越远的地方投影误差越大。
``` 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表示投影坐标，GEOGCS表示球面坐标。
解决方法，不要转换为26986坐标系统(即平面坐标)，如果求球面坐标，使用PostGIS对应函数即可。即使要求平面距离，那么也应该使用国内常用坐标系（使用央子午线在国内的坐标系，参考本文前面的内容）。
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)
```
PostgreSQL 许愿链接
您的愿望将传达给PG kernel hacker、数据库厂商等, 帮助提高数据库产品质量和功能, 说不定下一个PG版本就有您提出的功能点. 针对非常好的提议，奖励限量版PG文化衫、纪念品、贴纸、PG热门书籍等，奖品丰富，快来许愿。开不开森.
9.9元购买3个月阿里云RDS PostgreSQL实例
PostgreSQL 解决方案集合
德哥 / digoal's github - 公益是一辈子的事.
