目錄
需求和背景
在已有大量經緯度坐標點的情況下，給定一組經緯度如何快速定位到附近的POI有哪些？
現在使用經緯度轉geohash的算法，將二維的距離運算轉換為like前綴匹配。通過比較9位到5位前綴，來獲取附近5米到3km之內的坐標，為了尋求更快的定位方法，測試一下postgres的空間類型。
安裝插件postgis
先安裝了pg-10, 並且是通過yum安裝的。導入過repo.
檢查插件
yum search postgis
postgis-docs.x86_64 : Extra documentation for PostGIS
postgis-jdbc.x86_64 : The JDBC driver for PostGIS
postgis-utils.x86_64 : The utils for PostGIS
postgis23_10-client.x86_64 : Client tools and their libraries of PostGIS
postgis23_10-devel.x86_64 : Development headers and libraries for PostGIS
postgis23_10-docs.x86_64 : Extra documentation for PostGIS
postgis23_10-utils.x86_64 : The utils for PostGIS
postgis24_10-client.x86_64 : Client tools and their libraries of PostGIS
postgis24_10-debuginfo.x86_64 : Debug information for package postgis24_10
postgis24_10-devel.x86_64 : Development headers and libraries for PostGIS
postgis24_10-docs.x86_64 : Extra documentation for PostGIS
postgis24_10-utils.x86_64 : The utils for PostGIS
postgis.x86_64 : Geographic Information Systems Extensions to PostgreSQL
postgis23_10.x86_64 : Geographic Information Systems Extensions to PostgreSQL
postgis24_10.x86_64 : Geographic Information Systems Extensions to PostgreSQL
安裝
yum install postgis.x86_64 postgis24_10.x86_64
系統安裝了插件之后，數據庫還要繼續啟用插件才行。
針對數據庫啟用插件
# 添加空間插件
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;
安裝之后，public下會新增一個表spatial_ref_sys。
點POINT類型和距離
點POINT類型的數據結構為POINT(0 0)，正好可以用作存儲經緯度。
表添加POINT類型
AddGeometryColumn
使用函數AddGeometryColumn, 命令行查看函數
\df+ AddGeometryColumn
Synopsis
text AddGeometryColumn(varchar table_name, varchar column_name, integer srid, varchar type, integer dimension, boolean use_typmod=true);
text AddGeometryColumn(varchar schema_name, varchar table_name, varchar column_name, integer srid, varchar type, integer dimension, boolean use_typmod=true);
text AddGeometryColumn(varchar catalog_name, varchar schema_name, varchar table_name, varchar column_name, integer srid, varchar type, integer dimension, boolean use_typmod=true);
添加兩個點字段
SELECT AddGeometryColumn ('poi', 'geom_point', 4326, 'POINT', 2);
SELECT AddGeometryColumn ('poi', 'geom_point_26986', 26986, 'POINT', 2);
其中兩個重要的坐標體系
- 4326 \ GCS_WGS_1984 \ World Geodetic System (WGS)
- 26986 \ 美國馬薩諸塞州地方坐標系（區域坐標系）\ 投影坐標, 平面坐標
直接添加
ALTER TABLE poi ADD COLUMN geom_p_alter geometry(POINT,4326);
添加空間索引
CREATE INDEX idx_point
ON poi
USING gist(geom_point);
插入點
使用函數將文本轉換為幾何類型: ST_GeomFromText
sdx=# SELECT ST_GeomFromText('POINT(120.377041 36.066019)', 4326);
st_geomfromtext
----------------------------------------------------
0101000020E61000001310937021185E4012F5824F73084240
(1 row)
使用坐標轉換函數轉換坐標體系：ST_Transform
sdx=# SELECT ST_Transform(ST_GeomFromText('POINT(120.377041 36.066019)', 4326),26986)
st_transform
----------------------------------------------------
01010000206A690000B6A9B046D9615AC162C3613707DD6441
使用函數將幾何類型轉換為文本描述：ST_AsText
SELECT ST_AsText(ST_GeomFromText('POINT(120.377041 36.066019)', 4326));
st_astext
-----------------------------
POINT(120.377041 36.066019)
(1 row)
插入三個點
update poi set
geom_point=ST_GeomFromText('POINT(121.248642 31.380415)', 4326),
geom_point_26986=ST_Transform(ST_GeomFromText('POINT(121.248642 31.380415)', 4326),26986),
geom_p_alter=ST_GeomFromText('POINT(121.248642 31.380415)', 4326)
WHERE uuid='462745f185a349bbb8454f70d085baae';
SELECT geom_point,geom_point_26986,geom_p_alter from poi WHERE uuid='462745f185a349bbb8454f70d085baae';
geom_point | geom_point_26986 | geom_p_alter
----------------------------------------------------+----------------------------------------------------+--
0101000020E610000085766FC1E94F5E400D1AFA27B858D83F | 01010000206A69000087930146005A5CC1ECE89370F91A6541 | 0101000020E610000085766FC1E94F5E400D1AFA27B858D83F
(1 row)
驗證：
SELECT ST_AsText(geom_point),ST_AsText(geom_point_26986),geom_p_alter from s_poi_gaode WHERE uuid='462745f185a349bbb8454f70d085baae';
st_astext | st_astext | geom_p_alter
-----------------------------+-------------------------------------------+----------------------------------
POINT(121.248642 31.380415) | POINT(-7432193.09384621 11065291.5180554) | 0101000020E610000085766FC1E94F5E400D1AFA27B858D83F
(1 row)
批量更新現有的經緯度字段為POINT
update s_poi_gaode set
geom_point=ST_GeomFromText('POINT('||longitude||' ' ||latitude||')', 4326),
geom_point_26986=ST_Transform(ST_GeomFromText('POINT('||longitude||' ' ||latitude||')', 4326),26986);
驗證：
SELECT longitude,latitude,ST_AsText(geom_point),ST_AsText(geom_point_26986) from s_poi_gaode WHERE uuid='e3ebbcf15cc545408ac8b22d4df64ca6';
longitude | latitude | st_astext | st_astext
------------+-----------+-----------------------------+------------------------------------------
121.417666 | 31.281433 | POINT(121.417666 31.281433) | POINT(-7448729.03389232 11054385.435284)
(1 row)
其中，需要注意的是，使用pg的字符串拼接符號||，POINT經緯度之間要留空格。
兩個點之間的距離
距離計算函數
ST_Distance
文本轉換地理幾何類型函數
ST_GeogFromText 。
文本轉換為地理幾何類型函數
ST_GeographyFromText
計算距離，單位是m的方法
-- 921.37629155
select ST_Distance(ST_GeographyFromText('SRID=4326;POINT(114.017299 22.537126)'),
ST_GeographyFromText('SRID=4326;POINT(114.025919 22.534866)')
);
-- 921.37629155
SELECT ST_Distance(ST_GeomFromText('POINT(114.017299 22.537126)',4326):: geography,
ST_GeomFromText('POINT(114.025919 22.534866)', 4326):: geography
);
-- 920.28519
SELECT ST_DistanceSphere(ST_GeomFromText('POINT(114.017299 22.537126)',4326),
ST_GeomFromText('POINT(114.025919 22.534866)', 4326)
);
-- unit=m 26986 馬薩諸塞州 投影平面坐標系 單位m result=972.989337453172
SELECT ST_Distance(
ST_Transform(ST_GeomFromText('POINT(114.017299 22.537126)',4326 ),26986),
ST_Transform(ST_GeomFromText('POINT(114.025919 22.534866)', 4326 ),26986)
);
計算距離，單位是度
# unit=degrees result=0.00891134108875483
SELECT ST_Distance(ST_GeomFromText('POINT(114.017299 22.537126)',4326),
ST_GeomFromText('POINT(114.025919 22.534866)', 4326)
);
關於單位是m的, 前三種的計算結果是正確的。最后一種坐標轉換的計算方法，
參考PostGIS 坐標轉換(SRID)的邊界問題引發的專業知識 - ST_Transform 建議國內不要使用馬薩諸塞州的投影平面，會使得距離計算不夠准確。
附近5公里內的點
使用函數ST_DWithin 可以計算兩個點之間的距離是否在5公里內。
# 計算兩個點是否在給定距離內
# 單位米m
SELECT ST_DWithin(
ST_GeographyFromText('SRID=4326;POINT(114.017299 22.537126)'),
ST_GeographyFromText('SRID=4326;POINT(114.025919 22.534866)'),
1000);
# 單位度degrees
SELECT ST_DWithin(
ST_GeomFromText('POINT(114.017299 22.537126)',4326),
ST_GeomFromText('POINT(114.025919 22.534866)', 4326),
0.00811134108875483);
-- 查找給定經緯度5km以內的點
SELECT
uuid,
longitude,
latitude,
ST_DistanceSphere (
geom_point,
ST_GeomFromText ( 'POINT(121.248642 31.380415)', 4326 )) distance
FROM
s_poi_gaode
WHERE
ST_DWithin ( geom_point :: geography, ST_GeomFromText ( 'POINT(121.248642 31.380415)', 4326 ) :: geography, 5000 ) IS TRUE
order by distance desc
LIMIT 30;
通過指定類型geom_point :: geography，單位變成米, 否則默認距離單位是度。
最近的10個點
SELECT * FROM s_poi_gaode_gps ORDER BY geom_point <-> ST_GeomFromText ( 'POINT(121.248642 31.380415)', 4326 ) LIMIT 10;
速度極快。
面多邊形'POLYGON'
添加字段類型
SELECT AddGeometryColumn ('basic_mall_v1', 'geom_fence', 4326, 'POLYGON', 2);
或者
ALTER TABLE basic_mall_v1 ADD COLUMN geom_fence_alter geometry(POLYGON,4326);
添加索引
CREATE INDEX idx_area_fence
ON basic_mall_v1
USING gist(geom_fence);
插入值
使用函數 ST_GeomFromText
SELECT ST_GeomFromText ( 'POLYGON((118.902957 32.085437,118.9041 32.086069,118.904754 32.085219,118.903592 32.084564,118.902957 32.085437))', 4326 );
st_geomfromtext
-----------------------------------------------------------------------
0103000020E610000001000000050000006F2C280CCAB95D40266F8099EF0A404012143FC6DCB95D4087191A4F040B4040363B527DE7B95D40B9FFC874E80A4040583B8A73D4B95D40A0353FFED20A40406F2C280CCAB95D40266F8099EF0A4040
(1 row)
-- 驗證
SELECT
ST_AsText(
ST_GeomFromText ( 'POLYGON((118.902957 32.085437,118.9041 32.086069,118.904754 32.085219,118.903592 32.084564,118.902957 32.085437))', 4326 )
)
st_astext
-------------------------------------------------------------------------------------------------------------------
POLYGON((118.902957 32.085437,118.9041 32.086069,118.904754 32.085219,118.903592 32.084564,118.902957 32.085437))
(1 row)
更新到數據庫字段
UPDATE basic_mall_v1 SET geom_fence=ST_GeomFromText ( 'POLYGON((118.902957 32.085437,118.9041 32.086069,118.904754 32.085219,118.903592 32.084564,118.902957 32.085437))', 4326 ) WHERE id=1000001;
-- 驗證：
SELECT ST_AsText(geom_fence) FROM basic_mall_v1 WHERE id=1000001;
st_astext
-------------------------------------------------------------------------------------------------------------------
POLYGON((118.902957 32.085437,118.9041 32.086069,118.904754 32.085219,118.903592 32.084564,118.902957 32.085437))
(1 row)
實際上，我們原始圍欄數據可能是這樣的
-- longitude,latitude; longitude,latitude; longitude,latitude;
119.306413,26.131464;119.307575,26.131739;119.30776,26.131224;119.307336,26.131114;119.307438,26.130791;119.306776,26.13059;119.306413,26.131464
需要將這個字段轉換成空間類型的圍欄字段。
UPDATE basic_mall_v1 SET geom_fence=
ST_GeomFromText(
'POLYGON(('||
replace(
replace(gaode_shape, ',', ' ' ),
';', ',')
||'))'
, 4326)
計算gps附近30m內的圍欄
使用函數ST_DWithin 判斷一個幾何對象是否在另一個的r距離以內：
SELECT
ST_Distance(ST_GeomFromText('POINT(120.731069 30.758984)',4326):: geography,
geom_fence :: geography
) AS distance, id, name
FROM
basic_mall_v1
WHERE
ST_DWithin (
geom_fence :: geography,
ST_GeomFromText ( 'POINT( 120.731069 30.758984)', 4326 ) :: geography,
30
) ORDER BY distance
LIMIT 10;
使用函數boolean ST_Within(geometry A, geometry B); 判斷A是否完全在B內部
SELECT
id, name
FROM
basic_mall_v1
WHERE
ST_Within (
ST_GeomFromText('POINT('|| #{longitude} ||' '|| #{latitude} ||')',4326),
geom_fence
)
更多測試見菜鳥末端軌跡 - 電子圍欄(解密支撐每天251億個包裹的數據庫)
關於坐標體系
地理坐標系（Geographic coordinate system）
首先理解地理坐標系（Geographic coordinate system），Geographic coordinate system直譯為地理坐標系統，是以經緯度為地圖的存儲單位的。
很明顯，Geographic coordinate system是球面坐標系統。我們要將地球上的數字化信息存放到球面坐標系統上，如何進行操作呢？地球是一個不規則的橢球，如何將數據信息以科學的方法存放到橢球上？
這必然要求我們找到這樣的一個橢球體。這樣的橢球體具有特點：
可以量化計算的。具有長半軸，短半軸，偏心率。
以下幾行便是Krasovsky_1940橢球及其相應參數。
Alias:
Abbreviation:
Remarks:
Angular Unit: Degree (0.017453292519943299)
Prime Meridian（起始經度）: Greenwich (0.000000000000000000)
Datum（大地基准面）: D_Beijing_1954
Spheroid（參考橢球體）: Krasovsky_1940
Semimajor Axis: 6378245.000000000000000000
Semiminor Axis: 6356863.018773047300000000
Inverse Flattening: 298.300000000000010000
然而有了這個橢球體以后還不夠，還需要一個大地基准面將這個橢球定位。在坐標系統描述中，可以看到有這么一行：
Datum: D_Beijing_1954
表示，大地基准面是D_Beijing_1954。
Projection coordinate system（投影坐標系統）
投影坐標系統，實質上便是平面坐標系統，其地圖單位通常為米.
投影的意義：將球面坐標轉化為平面坐標的過程便稱為投影。
參數
Projection: Gauss_Kruger
Parameters:
False_Easting: 500000.000000
False_Northing: 0.000000
Central_Meridian: 117.000000
Scale_Factor: 1.000000
Latitude_Of_Origin: 0.000000
Linear Unit: Meter (1.000000)
Geographic Coordinate System:
Name: GCS_Beijing_1954
Alias:
Abbreviation:
Remarks:
Angular Unit: Degree (0.017453292519943299)
Prime Meridian: Greenwich (0.000000000000000000)
Datum: D_Beijing_1954
Spheroid: Krasovsky_1940
Semimajor Axis: 6378245.000000000000000000
Semiminor Axis: 6356863.018773047300000000
Inverse Flattening: 298.300000000000010000