102
PostGIS 距離計算規範 - 投影 與 球 坐標係, geometry 與 geography 類型
標簽
PostgreSQL , PostGIS , 球坐標 , 平麵坐標 , 球麵距離 , 平麵距離
背景
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、計算球麵距離
https://postgis.net/docs/manual-2.4/ST_DistanceSphere.html
2、計算球麵以及平麵距離(取決於輸入的類型geometry還是geography)
https://postgis.net/docs/manual-2.4/ST_Distance.html
3、坐標係轉換
https://postgis.net/docs/manual-2.4/ST_Transform.html
4、投影坐標和球坐標
5、PostGIS學習建議
6、https://blog.163.com/jey_df/blog/static/18255016120149145755781/
https://desktop.arcgis.com/zh-cn/arcmap/10.3/guide-books/map-projections/mercator.htm
最後更新：2017-10-28 23:34:38