PostGIS 距离计算建议 - 投影 与 球 坐标系, geometry 与 geography 类型
背景
PostGIS中有两种常用的空间类型geometry和geography，这两种数据类型有什么差异，应该如何选择？
对于GIS来说，首先是坐标系，有两种：一种是球坐标（地理坐标），另一种是平面坐标（投影坐标）。
球坐标通常用于计算，平面坐标通常用于展示（也可以计算）。
投影坐标是从球坐标投影后展开得来(用一个圆柱将地球包起来，把地球当成会发光的光源，投影后，将圆柱展开得到)，投影的范围越大，精度就越低。范围越小，精度就越高。除了投影扇形的大小有区别，在不同的行业，也有不同的坐标系，例如用于测绘的坐标系。
目前用得最多的有SRID=4326球坐标，SRID为EPSG:3785的墨卡托投影坐标。
再来说geometry和geography两种类型，geometry支持平面对象也支持空间对象，而geography则仅支持空间对象。
geometry支持更多的函数，一些几何计算的代价更低。
geography支持的函数略少，计算代价更高。那为什么还要geography呢？因
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”
既然提到距离计算和投影坐标系有关，引入了本文的问题：
在不知道要计算的geometry点，在什么投影坐标系下时，往往计算距离得到的结果并不准确。
例如，下面的点，换算成2163坐标系，计算距离得到的结果并不准确。
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)
2163坐标系内容如下：
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
这样会造成一个假象如下：
用st_distance函数计算出来的经纬度之间的距离，跟用java程序算出来的距离相差很大。
这个例子，st_distance算出来的距离是4030，我们程序算出来的是4445，换另外两个相距很远的点，这个差距值会更大。
正确姿势
算球面距离，不要算直线距离。
1、使用球坐标，通过st_distancespheroid计算两个geometry类型的球面距离。
对于geometry类型，可以使用st_distancespheroid直接用球坐标计算，在计算时会自动设置这个椭球特性（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、使用平面坐标，通过st_distance计算两个geometry类型的平面投影后的距离。
采用精准的投影坐标（小面积投影坐标系）（但是必须要覆盖到要计算的两个点）
ST_Transform(ST_GeomFromText('POINT(120.08 30.96)', 4326), 2163 )
如果geometry值的SRID不是（高精度）目标坐标系，可以使用ST_Transform函数进行转换，转换为目标投影坐标系，再计算距离。
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、使用球坐标，通过ST_Distance计算两个geography类型的球面距离。
float ST_Distance(geography gg1, geography gg2, boolean use_spheroid); -- 适用椭球体（WGS84）
use_spheroid设置为ture表示使用:
-- WGS84 椭球体参数定义
vspheroid := 'SPHEROID["WGS84",6378137,298.257223563]' ;
这里的XXXX就是你要选择的球坐标系SRID。在spatial_ref_sys表里可以查看各种坐标系。
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)
如果允许一定的偏差，可以使用全球球坐标系4326。
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)
geography只支持球坐标系，使用投影坐标会报错。
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、指定SPHEROID内容，通过st_distancesphereoid计算geometry的球面距离。
这种方法最为精确，但是要求了解计算距离当地的地形属性（即输入参数的spheroid的内容）。
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)
小结
计算距离，应该考虑到被计算的两点所在处的地球特性（spheroid）。这样计算得到的距离才是最精确的。
geometry和geography类型的选择，建议使用geometry，既能支持球坐标系，又能支持平面坐标系。主要考虑到用户是否了解位置所在处的地理特性，选择合适的坐标系。
-- 适用平面直角坐标，适用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);
参考
1、计算球面距离
http://postgis.net/docs/manual-2.4/ST_DistanceSphere.html
2、计算球面以及平面距离(取决于输入的类型geometry还是geography)
http://postgis.net/docs/manual-2.4/ST_Distance.html
3、坐标系转换
http://postgis.net/docs/manual-2.4/ST_Transform.html
4、投影坐标和球坐标
5、PostGIS学习建议
6、http://blog.163.com/jey_df/blog/static/18255016120149145755781/
https://desktop.arcgis.com/zh-cn/arcmap/10.3/guide-books/map-projections/mercator.htm