PostGIS 距离计算规范 - 投影与球坐标系, geometry 与 geography 类型

清疚 2022-10-12 01:02 445阅读 0赞

平面坐标，球面坐标

SRID=2369  投影坐标  
SRID=4326  球面坐标  
ST\_Transform(ST\_GeomFromText('POINT(120.08 30.96)', 4326), 2390 )  坐标系转换函数

PostGIS中有两种常用的空间类型geometry和geography，这两种数据类型有什么差异，应该如何选择？

对于GIS来说，首先是坐标系，有两种：一种是球坐标（地理坐标），另一种是平面坐标（投影坐标）。

球坐标通常用于计算，平面坐标通常用于展示（也可以计算）。

投影坐标是从球坐标投影后展开得来(用一个圆柱将地球包起来，把地球当成会发光的光源，投影后，将圆柱展开得到)，投影的范围越大，精度就越低。范围越小，精度就越高。除了投影扇形的大小有区别，在不同的行业，也有不同的坐标系，例如用于测绘的坐标系。

目前用得最多的有SRID=4326球坐标，SRID为EPSG:3785的墨卡托投影坐标。

再来说geometry和geography两种类型，geometry支持平面对象也支持空间对象，而geography则仅支持空间对象。

geometry支持更多的函数，一些几何计算的代价更低。

geography支持的函数略少，计算代价更高。那为什么还要geography呢？因

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

既然提到距离计算和投影坐标系有关，引入了本文的问题：

在不知道要计算的geometry点，在什么投影坐标系下时，往往计算距离得到的结果并不准确。

例如，下面的点，换算成2163坐标系，计算距离得到的结果并不准确。

1.  `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 ));`
2.  `st_distance`
3.  `------------------`
4.  `4030.46766234184`
5.  `(1 row)`

2163坐标系内容如下：

1.  `postgres=# select * from spatial_ref_sys where srid=2163;`
2.  `-[ RECORD 1 ]-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------`
3.  `srid | 2163`
4.  `auth_name | EPSG`
5.  `auth_srid | 2163`
6.  `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"]]`
7.  `proj4text | +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs`
8.

1.  `AUTHORITY["EPSG", "9122"]指的是EPSG数据集中UNIT为degree的ID是9122；`
2.  
3.  `AUTHORITY["EPSG", "4326"]指的是地理坐标系WGS 84的ID是4326；`
4.  
5.  `AUTHORITY["EPSG", "9001"]指的是EPSG中UNIT为meter的ID是9001；`
6.  
7.  `AUTHORITY["EPSG", "32650"]指的是该投影坐标系WGS 84 / UTM zone 50N的ID是32650`

这样会造成一个假象如下：

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

## 正确姿势 ##

算球面距离，不要算直线距离。

### 1、使用球坐标，通过st\_distancespheroid计算两个geometry类型的球面距离。 ###

对于geometry类型，可以使用st\_distancespheroid直接用球坐标计算，在计算时会自动设置这个椭球特性（SPHEROID\["Krasovsky\_1940",6378245.000000,298.299997264589\] ）。

1.  `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]');`
2.  `st_distancespheroid`
3.  `---------------------`
4.  `4434.73734584354`
5.  `(1 row)`

### 2、使用平面坐标，通过st\_distance计算两个geometry类型的平面投影后的距离。 ###

采用精准的投影坐标（小面积投影坐标系）（但是必须要覆盖到要计算的两个点）

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

如果geometry值的SRID不是（高精度）目标坐标系，可以使用ST\_Transform函数进行转换，转换为目标投影坐标系，再计算距离。

1.  `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 ));`
2.  `st_distance`
3.  `------------------`
4.  `4547.55477647394`
5.  `(1 row)`

### 3、使用球坐标，通过ST\_Distance计算两个geography类型的球面距离。 ###

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

这里的XXXX就是你要选择的球坐标系SRID。在spatial\_ref\_sys表里可以查看各种坐标系。

1.  `postgres=# SELECT st_distance(ST_GeogFromText('SRID=xxxx;POINT(120.08 30.96)'), ST_GeogFromText('SRID=xxxx;POINT(120.08 30.92)'), true);`
2.  `st_distance`
3.  `----------------`
4.  `xxxxxxxxxxxxxxxxxxxx`
5.  `(1 row)`
6.  
7.  `例如`
8.  
9.  `postgres=# SELECT st_distance(ST_GeogFromText('SRID=4610;POINT(120.08 30.96)'), ST_GeogFromText('SRID=4610;POINT(120.08 30.92)'), true);`
10. `st_distance`
11. `----------------`
12. `4434.739418211`
13. `(1 row)`

如果允许一定的偏差，可以使用全球球坐标系4326。

1.  `postgres=# SELECT st_distance(ST_GeogFromText('SRID=4326;POINT(120.08 30.96)'), ST_GeogFromText('SRID=4326;POINT(120.08 30.92)'), true);`
2.  `st_distance`
3.  `----------------`
4.  `4434.737345844`
5.  `(1 row)`

geography只支持球坐标系，使用投影坐标会报错。

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

### 4、指定SPHEROID内容，通过st\_distancesphereoid计算geometry的球面距离。 ###

这种方法最为精确，但是要求了解计算距离当地的地形属性（即输入参数的spheroid的内容）。

1.  `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]');`
2.  `st_distancesphere`
3.  `-------------------`
4.  `4447.803189385`
5.  `(1 row)`

## 小结 ##

计算距离，应该考虑到被计算的两点所在处的地球特性（spheroid）。这样计算得到的距离才是最精确的。

geometry和geography类型的选择，建议使用geometry，既能支持球坐标系，又能支持平面坐标系。主要考虑到用户是否了解位置所在处的地理特性，选择合适的坐标系。

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

## 参考 ##

1、计算球面距离

[http://postgis.net/docs/manual-2.4/ST\_DistanceSphere.html][http_postgis.net_docs_manual-2.4_ST_DistanceSphere.html]

2、计算球面以及平面距离(取决于输入的类型geometry还是geography)

[http://postgis.net/docs/manual-2.4/ST\_Distance.html][http_postgis.net_docs_manual-2.4_ST_Distance.html]

3、坐标系转换

[http://postgis.net/docs/manual-2.4/ST\_Transform.html][http_postgis.net_docs_manual-2.4_ST_Transform.html]

4、投影坐标和球坐标

[《地理坐标系（球面坐标系）和投影坐标系（平面坐标系）》][Link 1]

5、PostGIS学习建议

[《PostGIS 空间数据学习建议》][PostGIS]

6、[http://blog.163.com/jey\_df/blog/static/18255016120149145755781/][http_blog.163.com_jey_df_blog_static_18255016120149145755781]

[https://desktop.arcgis.com/zh-cn/arcmap/10.3/guide-books/map-projections/mercator.htm][https_desktop.arcgis.com_zh-cn_arcmap_10.3_guide-books_map-projections_mercator.htm]

[http_postgis.net_docs_manual-2.4_ST_DistanceSphere.html]: http://postgis.net/docs/manual-2.4/ST_DistanceSphere.html
[http_postgis.net_docs_manual-2.4_ST_Distance.html]: http://postgis.net/docs/manual-2.4/ST_Distance.html
[http_postgis.net_docs_manual-2.4_ST_Transform.html]: http://postgis.net/docs/manual-2.4/ST_Transform.html
[Link 1]: https://github.com/digoal/blog/blob/master/201709/20170911_01.md
[PostGIS]: https://github.com/digoal/blog/blob/master/201708/20170809_01.md
[http_blog.163.com_jey_df_blog_static_18255016120149145755781]: http://blog.163.com/jey_df/blog/static/18255016120149145755781/
[https_desktop.arcgis.com_zh-cn_arcmap_10.3_guide-books_map-projections_mercator.htm]: https://desktop.arcgis.com/zh-cn/arcmap/10.3/guide-books/map-projections/mercator.htm