1/* contrib/earthdistance/earthdistance--1.1.sql */ 2 3-- complain if script is sourced in psql, rather than via CREATE EXTENSION 4\echo Use "CREATE EXTENSION earthdistance" to load this file. \quit 5 6-- earth() returns the radius of the earth in meters. This is the only 7-- place you need to change things for the cube base distance functions 8-- in order to use different units (or a better value for the Earth's radius). 9 10CREATE FUNCTION earth() RETURNS float8 11LANGUAGE SQL IMMUTABLE PARALLEL SAFE 12AS 'SELECT ''6378168''::float8'; 13 14-- Astronomers may want to change the earth function so that distances will be 15-- returned in degrees. To do this comment out the above definition and 16-- uncomment the one below. Note that doing this will break the regression 17-- tests. 18-- 19-- CREATE FUNCTION earth() RETURNS float8 20-- LANGUAGE SQL IMMUTABLE 21-- AS 'SELECT 180/pi()'; 22 23-- Define domain for locations on the surface of the earth using a cube 24-- datatype with constraints. cube provides 3D indexing. 25-- The cube is restricted to be a point, no more than 3 dimensions 26-- (for less than 3 dimensions 0 is assumed for the missing coordinates) 27-- and that the point must be very near the surface of the sphere 28-- centered about the origin with the radius of the earth. 29 30CREATE DOMAIN earth AS cube 31 CONSTRAINT not_point check(cube_is_point(value)) 32 CONSTRAINT not_3d check(cube_dim(value) <= 3) 33 CONSTRAINT on_surface check(abs(cube_distance(value, '(0)'::cube) / 34 earth() - '1'::float8) < '10e-7'::float8); 35 36CREATE FUNCTION sec_to_gc(float8) 37RETURNS float8 38LANGUAGE SQL 39IMMUTABLE STRICT 40PARALLEL SAFE 41AS 'SELECT CASE WHEN $1 < 0 THEN 0::float8 WHEN $1/(2*earth()) > 1 THEN pi()*earth() ELSE 2*earth()*asin($1/(2*earth())) END'; 42 43CREATE FUNCTION gc_to_sec(float8) 44RETURNS float8 45LANGUAGE SQL 46IMMUTABLE STRICT 47PARALLEL SAFE 48AS 'SELECT CASE WHEN $1 < 0 THEN 0::float8 WHEN $1/earth() > pi() THEN 2*earth() ELSE 2*earth()*sin($1/(2*earth())) END'; 49 50CREATE FUNCTION ll_to_earth(float8, float8) 51RETURNS earth 52LANGUAGE SQL 53IMMUTABLE STRICT 54PARALLEL SAFE 55AS 'SELECT cube(cube(cube(earth()*cos(radians($1))*cos(radians($2))),earth()*cos(radians($1))*sin(radians($2))),earth()*sin(radians($1)))::earth'; 56 57CREATE FUNCTION latitude(earth) 58RETURNS float8 59LANGUAGE SQL 60IMMUTABLE STRICT 61PARALLEL SAFE 62AS 'SELECT CASE WHEN cube_ll_coord($1, 3)/earth() < -1 THEN -90::float8 WHEN cube_ll_coord($1, 3)/earth() > 1 THEN 90::float8 ELSE degrees(asin(cube_ll_coord($1, 3)/earth())) END'; 63 64CREATE FUNCTION longitude(earth) 65RETURNS float8 66LANGUAGE SQL 67IMMUTABLE STRICT 68PARALLEL SAFE 69AS 'SELECT degrees(atan2(cube_ll_coord($1, 2), cube_ll_coord($1, 1)))'; 70 71CREATE FUNCTION earth_distance(earth, earth) 72RETURNS float8 73LANGUAGE SQL 74IMMUTABLE STRICT 75PARALLEL SAFE 76AS 'SELECT sec_to_gc(cube_distance($1, $2))'; 77 78CREATE FUNCTION earth_box(earth, float8) 79RETURNS cube 80LANGUAGE SQL 81IMMUTABLE STRICT 82PARALLEL SAFE 83AS 'SELECT cube_enlarge($1, gc_to_sec($2), 3)'; 84 85--------------- geo_distance 86 87CREATE FUNCTION geo_distance (point, point) 88RETURNS float8 89LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE AS 'MODULE_PATHNAME'; 90 91--------------- geo_distance as operator <@> 92 93CREATE OPERATOR <@> ( 94 LEFTARG = point, 95 RIGHTARG = point, 96 PROCEDURE = geo_distance, 97 COMMUTATOR = <@> 98); 99