Postgres2015全国用户大会将于11月20至21日在北京丽亭华苑酒店召开。本次大会嘉宾阵容强大，国内顶级PostgreSQL数据库专家将悉数到场，并特邀欧洲、俄罗斯、日本、美国等国家和地区的数据库方面专家助阵:
- Postgres-XC项目的发起人铃木市一(SUZUKI Koichi)
- Postgres-XL的项目发起人Mason Sharp
- pgpool的作者石井达夫(Tatsuo Ishii)
- PG-Strom的作者海外浩平(Kaigai Kohei)
- Greenplum研发总监姚延栋
- 周正中(德哥), PostgreSQL中国用户会创始人之一
- 汪洋，平安科技数据库技术部经理
- ……
|
 |
PostgreSQL提供了一个扩展模块earthdistance，实际上是将地球构造为一个标准的圆球体（实际上是扁球体），利用cube或point来表示地球上的点。
其中cube是用来记录球坐标的，通过坐标来表示地球上的点。（实际上是做了一定约束的cube, 即域类型earth）
-- Define domain for locations on the surface of the earth using a cube-- datatype with constraints. cube provides 3D indexing.-- The cube is restricted to be a point, no more than 3 dimensions-- (for less than 3 dimensions 0 is assumed for the missing coordinates)-- and that the point must be very near the surface of the sphere-- centered about the origin with the radius of the earth.
有三个约束CREATE DOMAIN earth AS cubeCONSTRAINT not_point check(cube_is_point(value))CONSTRAINT not_3d check(cube_dim(value) <= 3)CONSTRAINT on_surface check(abs(cube_distance(value, '(0)'::cube) /earth() - 1) < '10e-7'::float8);
而point是用来记录纬度和经度的，通过纬度和经度来表示地球上的点。
使用纬度和经度来表示有一定的缺陷，例如有正负，有边界需要注意，用球坐标则不存在这个问题。
所以需要依赖另外一个模块cube，这个模块是多维数据结构模型，可以表示多维的点，多维的BOX区域。
例如我要计算北京和杭州的地表距离，只需要提供两个位置经纬度就可以计算，使用earthdistance计算的结果没有使用postgis计算的结果精确，（原因地球是扁球体，而非标准圆球体）PostGIS考虑了这一点。
例子：
create extension cube;create extension earthdistance;
earth()函数，返回地球半径，单位千米。
CREATE FUNCTION earth() RETURNS float8LANGUAGE SQL IMMUTABLEAS 'SELECT ''6378168''::float8';
创建完后会建立一个域类型earth, 其实是基于cube建的。
postgres=# select * from pg_type where typname='earth';-[ RECORD 1 ]--+------------typname | earthtypnamespace | 2200typowner | 10typlen | -1typbyval | ftyptype | dtypcategory | Utypispreferred | ftypisdefined | ttypdelim | ,typrelid | 0typelem | 0typarray | 0typinput | domain_intypoutput | cube_outtypreceive | domain_recvtypsend | -typmodin | -typmodout | -typanalyze | -typalign | dtypstorage | ptypnotnull | ftypbasetype | 16545typtypmod | -1typndims | 0typcollation | 0typdefaultbin |typdefault |typacl |
将纬度，经度转换为earth坐标类型。第一个参数表示纬度，第二个参数表示经度。
postgres=# select ll_to_earth(30.3,120.2);-[ RECORD 1 ]--------------------------------------------------------ll_to_earth | (-2770071.42546437, 4759459.23949754, 3217961.94533299)
将球坐标转换为经度或纬度。
postgres=# select latitude(ll_to_earth(30.3,120.2));-[ RECORD 1 ]--latitude | 30.3postgres=# select longitude(ll_to_earth(30.3,120.2));-[ RECORD 1 ]----longitude | 120.2
直线距离和球体表面距离的换算。
直线距离转换为球面距离，
postgres=# select sec_to_gc(100000);-[ RECORD 1 ]--------------sec_to_gc | 100001.02425681
球面距离转换为直线距离，
postgres=# select gc_to_sec(100001.02425681);-[ RECORD 1 ]-----gc_to_sec | 100000
使用公式换算，用到圆周率pi()和地球半径：
CREATE FUNCTION sec_to_gc(float8)RETURNS float8LANGUAGE SQLIMMUTABLE STRICTAS 'SELECT CASE WHEN $1 < 0 THEN 0::float8 WHEN $1/(2*earth()) > 1 THEN pi()*earth() ELSE 2*earth()*asin($1/(2*earth())) END';
CREATE FUNCTION gc_to_sec(float8)RETURNS float8LANGUAGE SQLIMMUTABLE STRICTAS 'SELECT CASE WHEN $1 < 0 THEN 0::float8 WHEN $1/earth() > pi() THEN 2*earth() ELSE 2*earth()*sin($1/(2*earth())) END';
计算两个坐标的球面距离：
例如北京和杭州的球面距离，约1123公里。
postgres=# select earth_distance(ll_to_earth(39.9,116.4),ll_to_earth(30.3,120.2));-[ RECORD 1 ]--+----------------earth_distance | 1122999.5704515
算法，用到了CUBE的
cube_distance函数，先计算直线距离，再转换为球体表面距离
：
CREATE FUNCTION earth_distance(earth, earth)RETURNS float8LANGUAGE SQLIMMUTABLE STRICTAS 'SELECT sec_to_gc(cube_distance($1, $2))';
最后一个例子是范围，用户提供球体表面的一个坐标，以及一个半径信息，用来表示球体表面的一个以这个坐标为中心辐射的半径范围，但实际测试发现有一定的偏差，可能和cube的计算方法有关系。
earth_box(earth, float8) | cube | Returns a box suitable for an indexed search using the cube @> operator for points within a given great circle distance of a location. Some points in this box are further than the specified great circle distance from the location, so a second check using earth_distance should be included in the query. |
算法，还是要用到CUBE提供的
cube_enlarge函数，这个函数是将多维坐标的每个坐标都延展一个半径的距离，实际上已经不是球体了，而是立方体
：
CREATE FUNCTION earth_box(earth, float8)RETURNS cubeLANGUAGE SQLIMMUTABLE STRICTAS 'SELECT cube_enlarge($1, gc_to_sec($2), 3)';
第一个参数是坐标，第二个参数是直线距离（球面距离转换为直线距离），第三个参数是第一个参数的维度，这里是3维。所以用3表示。
| cube_enlarge(cube c, double r, int n) returns cube | Increases the size of a cube by a specified radius in at least n dimensions. If the radius is negative the cube is shrunk instead. This is useful for creating bounding boxes around a point for searching for nearby points. All defined dimensions are changed by the radius r. LL coordinates are decreased by r and UR coordinates are increased by r. If a LL coordinate is increased to larger than the corresponding UR coordinate (this can only happen when r < 0) than both coordinates are set to their average. If n is greater than the number of defined dimensions and the cube is being increased (r >= 0) then 0 is used as the base for the extra coordinates. |
例子：
将一个平面坐标(1,1)的每个方向扩展1个距离单位，得到
(0, 0),(2, 2)。你在纸上画一下就知道怎么回事了。
postgres=# select cube_enlarge('(1,1)',1,2);cube_enlarge---------------(0, 0),(2, 2)(1 row)
现在回到earthdistance的例子，例如我们前面计算得到杭州到北京的球面距离是1123公里，那么我以杭州为中心，辐射多少距离才能包含北京呢？如果是圆的话，辐射半径应该要达到
1123公里才能包含北京。但是实际上使用
earth_box得到的结果
并非如此。
先看看扩展10米得到的BOX，每个坐标方向正反各延展了10米，得到2个点，即一个BOX：
postgres=# select ll_to_earth(30.3,120.2);ll_to_earth---------------------------------------------------------(-2770071.42546437, 4759459.23949754, 3217961.94533299)(1 row)postgres=# select earth_box(ll_to_earth(30.3,120.2), 10);earth_box-----------------------------------------------------------------------------------------------------------------(-2770081.42546437, 4759449.23949754, 3217951.94533299),(-2770061.42546437, 4759469.23949754, 3217971.94533299)(1 row)
postgres=# select earth_box(ll_to_earth(30.3,120.2), 1122999.5704515);earth_box-----------------------------------------------------------------------------------------------------------------(-3891620.99816939, 3637909.66679251, 2096412.37262797),(-1648521.85275934, 5881008.81220257, 4339511.51803802)(1 row)
这个BOX是否包含北京呢？
postgres=# select earth_box(ll_to_earth(30.3,120.2), 1122999.5704515) @> ll_to_earth(39.9,116.4);?column?----------t(1 row)postgres=# select ll_to_earth(39.9,116.4);ll_to_earth--------------------------------------------------------(-2175648.05107066, 4382814.5785906, 4091273.51368619)(1 row)
实际上，800多公里就包含北京了，因为这个BOX是立方体，而不是球体，北京被包含在800多公里的球体外面，立方体里面。
postgres=# select earth_box(ll_to_earth(30.3,120.2), 880999.5704515) @> ll_to_earth(39.9,116.4);?column?----------t(1 row)
用PostGIS可以很好的解决这个问题。
earthdistance还支持使用point来表示坐标，计算得到的是英里单位。
经度，纬度。
postgres=# select point(116.4,39.9) <@> point(120.2,30.3);?column?-----------------697.01393638328(1 row)
转换为公里。postgres=# select 697.01393638328*1.60931;?column?-----------------------1121.7114979609763368(1 row)
与使用cube坐标计算得到的结果有一定的偏差，误差来自定义的地球半径不一样，见代码。
postgres=# select 3958.747716*1.60931;?column?------------------6370.85228683596(1 row)
earth()函数得到的是 AS 'SELECT ''6378168''::float8';
调整一下，两种坐标表示方法得到的结果就一致了：
postgres=# CREATE or REPLACE FUNCTION earth() RETURNS float8LANGUAGE SQL IMMUTABLEAS 'SELECT ''6370852''::float8';CREATE FUNCTION
postgres=# select earth_distance(ll_to_earth(39.9,116.4),ll_to_earth(30.3,120.2));earth_distance------------------1121711.44745797(1 row)
代码如下：
--------------- geo_distance as operator <@>CREATE OPERATOR <@> (LEFTARG = point,RIGHTARG = point,PROCEDURE = geo_distance,COMMUTATOR = <@>);
--------------- geo_distanceCREATE FUNCTION geo_distance (point, point)RETURNS float8LANGUAGE C IMMUTABLE STRICT AS 'MODULE_PATHNAME';
#include "postgres.h"
#include <math.h>
#include "utils/geo_decls.h" /* for Point */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
PG_MODULE_MAGIC;
/* Earth's radius is in statute miles. 英里单位 */
static const double EARTH_RADIUS = 3958.747716;
static const double TWO_PI = 2.0 * M_PI;
/******************************************************** geo_distance_internal - distance between points** args:* a pair of points - for each point,* x-coordinate is longitude in degrees west of Greenwich* y-coordinate is latitude in degrees above equator** returns: double* distance between the points in miles on earth's surface******************************************************/
static doublegeo_distance_internal(Point *pt1, Point *pt2){double long1,lat1,long2,lat2;double longdiff;double sino;
/* convert degrees to radians */
long1 = degtorad(pt1->x);lat1 = degtorad(pt1->y);
long2 = degtorad(pt2->x);lat2 = degtorad(pt2->y);
/* compute difference in longitudes - want < 180 degrees */longdiff = fabs(long1 - long2);if (longdiff > M_PI)longdiff = TWO_PI - longdiff;
sino = sqrt(sin(fabs(lat1 - lat2) / 2.) * sin(fabs(lat1 - lat2) / 2.) +cos(lat1) * cos(lat2) * sin(longdiff / 2.) * sin(longdiff / 2.));if (sino > 1.)sino = 1.;
return 2. * EARTH_RADIUS * asin(sino);}
1858
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
