]> granicus.if.org Git - postgresql/commitdiff
The attached patch defines functions for getting distances between
authorBruce Momjian <bruce@momjian.us>
Fri, 8 Nov 2002 20:20:22 +0000 (20:20 +0000)
committerBruce Momjian <bruce@momjian.us>
Fri, 8 Nov 2002 20:20:22 +0000 (20:20 +0000)
points on the surface of the earth and locating points within a
specified distance using an index based on the contrib/cube package. The
new functions are all of language type sql. A couple of bugs in the old
earthdistance function based on the point datatype are fixed. A
regression test has been added for both sets of functions. The README
file has been updated to include documentation on the new stuff. There
are comments about how this package is also useful for Astronomers.

Bruno Wolff III

contrib/earthdistance/README.earthdistance
contrib/earthdistance/earthdistance.c
contrib/earthdistance/earthdistance.sql.in
contrib/earthdistance/expected/earthdistance.out
contrib/earthdistance/sql/earthdistance.sql

index 2ed5c2816595c30c726bfbcff48a83c56a5fb284..913179d489b2a084d17166fd88fbec373cb141a4 100644 (file)
----------------------------------------------------------------------
-I corrected a bug in the geo_distance code where two double constants
-were declared as int. I changed the distance function to use the
-haversine formula which is more accurate for small distances.
-I added a regression test to the package. I added a grant statement
-to give execute access for geo_distance to public.
+This contrib package contains two different approaches to calculating
+great circle distances on the surface of the Earth. The one described
+first depends on the contrib/cube package (which MUST be installed before
+earthdistance is installed). The second one is based on the point
+datatype using latitude and longitude for the coordinates. The install
+script makes the defined functions executable by anyone.
+
+Make sure contrib/cube has been installed.
+make
+make install
+make installcheck
+
+To use these functions in a particular database as a postgres superuser do:
+psql databasename < earthdistance.sql
 
+-------------------------------------------
+contrib/cube based Earth distance functions
 Bruno Wolff III
 September 2002
+
+A spherical model of the Earth is used.
+
+Data is stored in cubes that are points (both corners are the same) using 3
+coordinates representing the distance from the center of the Earth.
+
+The radius of the Earth is obtained from the earth() function. It is
+given in meters. But by changing this one function you can change it
+to use some other units or to use a different value of the radius
+that you feel is more appropiate.
+
+This package also has applications to astronomical databases as well.
+Astronomers will probably want to change earth() to return a radius of
+180/pi() so that distances are in degrees.
+
+Functions are provided to allow for input in latitude and longitude (in
+degrees), to allow for output of latitude and longitude, to calculate
+the great circle distance between two points and to easily specify a
+bounding box usable for index searches.
+
+The functions are all 'sql' functions. If you want to make these functions
+executable by other people you will also have to make the referenced
+cube functions executable. cube(text), cube_distance(cube,cube),
+cube_ll_coord(cube,int) and cube_enlarge(cube,float8,int) are used indirectly
+by the earth distance functions. is_point(cube) and cube_dim(cube) are used
+in suggested constraints for data in domain earth. cube_ur_coord(cube,int)
+is used in the regression tests and might be useful for looking at bounding
+box coordinates in user applications.
+
+A domain of type cube named earth is defined. Since check constraints
+are not supported for domains yet, this isn't as useful as it might be.
+However the checks that should be supplied to all data of type earth are:
+
+constraint not_point check(is_point(earth))
+constraint not_3d check(cube_dim(earth) <= 3)
+constraint on_surface check(abs(cube_distance(earth, '(0)'::cube) /
+           earth() - 1) < '10e-12'::float8);
+
+The following functions are provided:
+
+earth() - Returns the radius of the earth in meters.
+
+sec_to_gc(float8) - Converts the normal straight line (secant) distance between
+between two points on the surface of the Earth to the great circle distance
+between them.
+
+gc_to_sec(float8) - Converts the great circle distance between two points
+on the surface of the Earth to the normal straight line (secant) distance
+between them.
+
+ll_to_cube(float8, float8) - Returns the location of a point on the surface of
+the Earth given its latitude (argument 1) and longitude (argument 2) in degrees.
+
+latitude(earth) - Returns the latitude in degrees of a point on the surface
+of the Earth.
+
+longitude(earth) - Returns the longitude in degrees of a point on the surface
+of the Earth.
+
+earth_distance(earth, earth) - Returns the great circle distance between
+two points on the surface of the Earth.
+
+earth_box(earth, float8)  - 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
+made at the same time.
+
+One advantage of using cube representation over a point using latitude and
+longitude for coordinates, is that you don't have to worry about special
+conditions at +/- 180 degrees of longitude or near the poles.
+
+Below is the documentation for the Earth distance operator that works
+with the point data type.
+
+---------------------------------------------------------------------
+
+I corrected a bug in the geo_distance code where two double constants
+were declared as int. I also changed the distance function to use
+the haversine formula which is more accurate for small distances.
+Bruno Wolff
+September 2002
+
 ---------------------------------------------------------------------
+
 Date: Wed, 1 Apr 1998 15:19:32 -0600 (CST)
 From: Hal Snyder <hal@vailsys.com>
 To: vmehr@ctp.com
index f069edf13102fde6b629e0e01a3141c4254bd1dd..19c81a5783fc2842149bf4f65317a335bd0ba115 100644 (file)
@@ -50,7 +50,7 @@ geo_distance(Point *pt1, Point *pt2)
                                long2,
                                lat2;
        double          longdiff;
-        double         sino;
+       double          sino;
        double     *resultp = palloc(sizeof(double));
 
        /* convert degrees to radians */
index 70dc32abab64425222efbd8e9c28458efb7d1b45..de381774ad02fcc77ca2c4a62e9edb3e8bb87e70 100644 (file)
@@ -3,6 +3,77 @@ SET search_path = public;
 
 SET autocommit TO 'on';
 
+-- The earth functions rely on contrib/cube having been installed and loaded.
+
+-- earth() returns the radius of the earth in meters. This is the only
+-- place you need to change things for the cube base distance functions
+-- in order to use different units (or a better value for the Earth's radius).
+
+CREATE OR REPLACE FUNCTION earth() RETURNS float8
+LANGUAGE 'sql' IMMUTABLE
+AS 'SELECT \'6378168\'::float8';
+
+-- Astromers may want to change the earth function so that distances will be
+-- returned in degrees. To do this comment out the above definition and
+-- uncomment the one below. Note that doing this will break the regression
+-- tests.
+--
+-- CREATE OR REPLACE FUNCTION earth() RETURNS float8
+-- LANGUAGE 'sql' IMMUTABLE
+-- AS 'SELECT 180/pi()';
+
+-- Define domain for locations on the surface of the earth using a cube
+-- datatype with constraints. cube provides 3D indexing.
+-- Check constraints aren't currently supported.
+
+CREATE DOMAIN earth AS cube;
+--  CONSTRAINT not_point check(is_point(earth))
+--  CONSTRAINT not_3d check(cube_dim(earth) <= 3)
+--  CONSTRAINT on_surface check(abs(cube_distance(earth, '(0)'::cube) /
+--  earth() - 1) < '10e-12'::float8);
+
+CREATE OR REPLACE FUNCTION sec_to_gc(float8) 
+RETURNS float8
+LANGUAGE 'sql'
+IMMUTABLE STRICT
+AS '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 OR REPLACE FUNCTION gc_to_sec(float8)
+RETURNS float8
+LANGUAGE 'sql'
+IMMUTABLE STRICT
+AS 'SELECT CASE WHEN $1 < 0 THEN 0::float8 WHEN $1/earth() > pi() THEN 2*earth() ELSE 2*earth()*sin($1/(2*earth())) END';
+
+CREATE OR REPLACE FUNCTION ll_to_earth(float8, float8)
+RETURNS earth
+LANGUAGE 'sql'
+IMMUTABLE STRICT
+AS 'SELECT cube(\'(\'||earth()*cos(radians($1))*cos(radians($2))||\',\'||earth()*cos(radians($1))*sin(radians($2))||\',\'||earth()*sin(radians($1))||\')\')';
+
+CREATE OR REPLACE FUNCTION latitude(earth)
+RETURNS float8
+LANGUAGE 'sql'
+IMMUTABLE STRICT
+AS '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';
+
+CREATE OR REPLACE FUNCTION longitude(earth)
+RETURNS float8
+LANGUAGE 'sql'
+IMMUTABLE STRICT
+AS 'SELECT degrees(atan2(cube_ll_coord($1, 2), cube_ll_coord($1, 1)))';
+
+CREATE OR REPLACE FUNCTION earth_distance(earth, earth)
+RETURNS float8
+LANGUAGE 'sql'
+IMMUTABLE STRICT
+AS 'SELECT sec_to_gc(cube_distance($1, $2))';
+
+CREATE OR REPLACE FUNCTION earth_box(earth, float8)
+RETURNS cube
+LANGUAGE 'sql'
+IMMUTABLE STRICT
+AS 'SELECT cube_enlarge($1, gc_to_sec($2), 3)';
+  
 --------------- geo_distance
 
 CREATE OR REPLACE FUNCTION geo_distance (point, point)
index b867205479cc55a63d7b037d524d5aba9f3ea278..845f4dec6d42c66c4e3e07047af9869abc66852b 100644 (file)
 --
 --
 -- first, define the datatype.  Turn off echoing so that expected file
--- does not depend on contents of earthdistance.sql.
+-- does not depend on contents of earthdistance.sql or cube.sql.
 --
 \set ECHO none
+psql:../cube/cube.sql:12: NOTICE:  ProcedureCreate: type cube is not yet defined
+psql:../cube/cube.sql:17: NOTICE:  Argument type "cube" is only a shell
+--
+-- The radius of the Earth we are using.
+--
+select earth()::numeric(20,5);
+     earth     
+---------------
+ 6378168.00000
+(1 row)
+
+--
+-- Convert straight line distances to great circle distances.
+--
+select (pi()*earth())::numeric(20,5);
+    numeric     
+----------------
+ 20037605.73216
+(1 row)
+
+select sec_to_gc(0)::numeric(20,5);
+ sec_to_gc 
+-----------
+   0.00000
+(1 row)
+
+select sec_to_gc(2*earth())::numeric(20,5);
+   sec_to_gc    
+----------------
+ 20037605.73216
+(1 row)
+
+select sec_to_gc(10*earth())::numeric(20,5);
+   sec_to_gc    
+----------------
+ 20037605.73216
+(1 row)
+
+select sec_to_gc(-earth())::numeric(20,5);
+ sec_to_gc 
+-----------
+   0.00000
+(1 row)
+
+select sec_to_gc(1000)::numeric(20,5);
+ sec_to_gc  
+------------
+ 1000.00000
+(1 row)
+
+select sec_to_gc(10000)::numeric(20,5);
+  sec_to_gc  
+-------------
+ 10000.00102
+(1 row)
+
+select sec_to_gc(100000)::numeric(20,5);
+  sec_to_gc   
+--------------
+ 100001.02426
+(1 row)
+
+select sec_to_gc(1000000)::numeric(20,5);
+   sec_to_gc   
+---------------
+ 1001027.07131
+(1 row)
+
+--
+-- Convert great circle distances to straight line distances.
+--
+select gc_to_sec(0)::numeric(20,5);
+ gc_to_sec 
+-----------
+   0.00000
+(1 row)
+
+select gc_to_sec(sec_to_gc(2*earth()))::numeric(20,5);
+   gc_to_sec    
+----------------
+ 12756336.00000
+(1 row)
+
+select gc_to_sec(10*earth())::numeric(20,5);
+   gc_to_sec    
+----------------
+ 12756336.00000
+(1 row)
+
+select gc_to_sec(pi()*earth())::numeric(20,5);
+   gc_to_sec    
+----------------
+ 12756336.00000
+(1 row)
+
+select gc_to_sec(-1000)::numeric(20,5);
+ gc_to_sec 
+-----------
+   0.00000
+(1 row)
+
+select gc_to_sec(1000)::numeric(20,5);
+ gc_to_sec  
+------------
+ 1000.00000
+(1 row)
+
+select gc_to_sec(10000)::numeric(20,5);
+ gc_to_sec  
+------------
+ 9999.99898
+(1 row)
+
+select gc_to_sec(100000)::numeric(20,5);
+  gc_to_sec  
+-------------
+ 99998.97577
+(1 row)
+
+select gc_to_sec(1000000)::numeric(20,5);
+  gc_to_sec   
+--------------
+ 998976.08618
+(1 row)
+
+--
+-- Set coordinates using latitude and longitude.
+-- Extract each coordinate separately so we can round them.
+--
+select cube_ll_coord(ll_to_earth(0,0),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,0),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,0),3)::numeric(20,5);
+ cube_ll_coord | cube_ll_coord | cube_ll_coord 
+---------------+---------------+---------------
+ 6378168.00000 |       0.00000 |       0.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(360,360),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(360,360),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(360,360),3)::numeric(20,5);
+ cube_ll_coord | cube_ll_coord | cube_ll_coord 
+---------------+---------------+---------------
+ 6378168.00000 |       0.00000 |       0.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(180,180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(180,180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(180,180),3)::numeric(20,5);
+ cube_ll_coord | cube_ll_coord | cube_ll_coord 
+---------------+---------------+---------------
+ 6378168.00000 |       0.00000 |       0.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(180,360),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(180,360),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(180,360),3)::numeric(20,5);
+ cube_ll_coord  | cube_ll_coord | cube_ll_coord 
+----------------+---------------+---------------
+ -6378168.00000 |       0.00000 |       0.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(-180,-360),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-180,-360),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-180,-360),3)::numeric(20,5);
+ cube_ll_coord  | cube_ll_coord | cube_ll_coord 
+----------------+---------------+---------------
+ -6378168.00000 |       0.00000 |       0.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(0,180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,180),3)::numeric(20,5);
+ cube_ll_coord  | cube_ll_coord | cube_ll_coord 
+----------------+---------------+---------------
+ -6378168.00000 |       0.00000 |       0.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(0,-180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,-180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,-180),3)::numeric(20,5);
+ cube_ll_coord  | cube_ll_coord | cube_ll_coord 
+----------------+---------------+---------------
+ -6378168.00000 |       0.00000 |       0.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(90,0),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(90,0),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(90,0),3)::numeric(20,5);
+ cube_ll_coord | cube_ll_coord | cube_ll_coord 
+---------------+---------------+---------------
+       0.00000 |       0.00000 | 6378168.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(90,180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(90,180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(90,180),3)::numeric(20,5);
+ cube_ll_coord | cube_ll_coord | cube_ll_coord 
+---------------+---------------+---------------
+       0.00000 |       0.00000 | 6378168.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(-90,0),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-90,0),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-90,0),3)::numeric(20,5);
+ cube_ll_coord | cube_ll_coord | cube_ll_coord  
+---------------+---------------+----------------
+       0.00000 |       0.00000 | -6378168.00000
+(1 row)
+
+select cube_ll_coord(ll_to_earth(-90,180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-90,180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-90,180),3)::numeric(20,5);
+ cube_ll_coord | cube_ll_coord | cube_ll_coord  
+---------------+---------------+----------------
+       0.00000 |       0.00000 | -6378168.00000
+(1 row)
+
+--
+-- Test getting the latitude of a location.
+--
+select latitude(ll_to_earth(0,0))::numeric(20,10);
+   latitude   
+--------------
+ 0.0000000000
+(1 row)
+
+select latitude(ll_to_earth(45,0))::numeric(20,10);
+   latitude    
+---------------
+ 45.0000000000
+(1 row)
+
+select latitude(ll_to_earth(90,0))::numeric(20,10);
+   latitude    
+---------------
+ 90.0000000000
+(1 row)
+
+select latitude(ll_to_earth(-45,0))::numeric(20,10);
+    latitude    
+----------------
+ -45.0000000000
+(1 row)
+
+select latitude(ll_to_earth(-90,0))::numeric(20,10);
+    latitude    
+----------------
+ -90.0000000000
+(1 row)
+
+select latitude(ll_to_earth(0,90))::numeric(20,10);
+   latitude   
+--------------
+ 0.0000000000
+(1 row)
+
+select latitude(ll_to_earth(45,90))::numeric(20,10);
+   latitude    
+---------------
+ 45.0000000000
+(1 row)
+
+select latitude(ll_to_earth(90,90))::numeric(20,10);
+   latitude    
+---------------
+ 90.0000000000
+(1 row)
+
+select latitude(ll_to_earth(-45,90))::numeric(20,10);
+    latitude    
+----------------
+ -45.0000000000
+(1 row)
+
+select latitude(ll_to_earth(-90,90))::numeric(20,10);
+    latitude    
+----------------
+ -90.0000000000
+(1 row)
+
+select latitude(ll_to_earth(0,180))::numeric(20,10);
+   latitude   
+--------------
+ 0.0000000000
+(1 row)
+
+select latitude(ll_to_earth(45,180))::numeric(20,10);
+   latitude    
+---------------
+ 45.0000000000
+(1 row)
+
+select latitude(ll_to_earth(90,180))::numeric(20,10);
+   latitude    
+---------------
+ 90.0000000000
+(1 row)
+
+select latitude(ll_to_earth(-45,180))::numeric(20,10);
+    latitude    
+----------------
+ -45.0000000000
+(1 row)
+
+select latitude(ll_to_earth(-90,180))::numeric(20,10);
+    latitude    
+----------------
+ -90.0000000000
+(1 row)
+
+select latitude(ll_to_earth(0,-90))::numeric(20,10);
+   latitude   
+--------------
+ 0.0000000000
+(1 row)
+
+select latitude(ll_to_earth(45,-90))::numeric(20,10);
+   latitude    
+---------------
+ 45.0000000000
+(1 row)
+
+select latitude(ll_to_earth(90,-90))::numeric(20,10);
+   latitude    
+---------------
+ 90.0000000000
+(1 row)
+
+select latitude(ll_to_earth(-45,-90))::numeric(20,10);
+    latitude    
+----------------
+ -45.0000000000
+(1 row)
+
+select latitude(ll_to_earth(-90,-90))::numeric(20,10);
+    latitude    
+----------------
+ -90.0000000000
+(1 row)
+
+--
+-- Test getting the longitude of a location.
+--
+select longitude(ll_to_earth(0,0))::numeric(20,10);
+  longitude   
+--------------
+ 0.0000000000
+(1 row)
+
+select longitude(ll_to_earth(45,0))::numeric(20,10);
+  longitude   
+--------------
+ 0.0000000000
+(1 row)
+
+select longitude(ll_to_earth(90,0))::numeric(20,10);
+  longitude   
+--------------
+ 0.0000000000
+(1 row)
+
+select longitude(ll_to_earth(-45,0))::numeric(20,10);
+  longitude   
+--------------
+ 0.0000000000
+(1 row)
+
+select longitude(ll_to_earth(-90,0))::numeric(20,10);
+  longitude   
+--------------
+ 0.0000000000
+(1 row)
+
+select longitude(ll_to_earth(0,90))::numeric(20,10);
+   longitude   
+---------------
+ 90.0000000000
+(1 row)
+
+select longitude(ll_to_earth(45,90))::numeric(20,10);
+   longitude   
+---------------
+ 90.0000000000
+(1 row)
+
+select longitude(ll_to_earth(90,90))::numeric(20,10);
+   longitude   
+---------------
+ 90.0000000000
+(1 row)
+
+select longitude(ll_to_earth(-45,90))::numeric(20,10);
+   longitude   
+---------------
+ 90.0000000000
+(1 row)
+
+select longitude(ll_to_earth(-90,90))::numeric(20,10);
+   longitude   
+---------------
+ 90.0000000000
+(1 row)
+
+select longitude(ll_to_earth(0,180))::numeric(20,10);
+   longitude    
+----------------
+ 180.0000000000
+(1 row)
+
+select longitude(ll_to_earth(45,180))::numeric(20,10);
+   longitude    
+----------------
+ 180.0000000000
+(1 row)
+
+select longitude(ll_to_earth(90,180))::numeric(20,10);
+   longitude    
+----------------
+ 180.0000000000
+(1 row)
+
+select longitude(ll_to_earth(-45,180))::numeric(20,10);
+   longitude    
+----------------
+ 180.0000000000
+(1 row)
+
+select longitude(ll_to_earth(-90,180))::numeric(20,10);
+   longitude    
+----------------
+ 180.0000000000
+(1 row)
+
+select longitude(ll_to_earth(0,-90))::numeric(20,10);
+   longitude    
+----------------
+ -90.0000000000
+(1 row)
+
+select longitude(ll_to_earth(45,-90))::numeric(20,10);
+   longitude    
+----------------
+ -90.0000000000
+(1 row)
+
+select longitude(ll_to_earth(90,-90))::numeric(20,10);
+   longitude    
+----------------
+ -90.0000000000
+(1 row)
+
+select longitude(ll_to_earth(-45,-90))::numeric(20,10);
+   longitude    
+----------------
+ -90.0000000000
+(1 row)
+
+select longitude(ll_to_earth(-90,-90))::numeric(20,10);
+   longitude    
+----------------
+ -90.0000000000
+(1 row)
+
+--
+-- For the distance tests the following is some real life data.
+--
+-- Chicago has a latitude of 41.8 and a longitude of 87.6.
+-- Albuquerque has a latitude of 35.1 and a longitude of 106.7.
+-- (Note that latitude and longitude are specified differently
+-- in the cube based functions than for the point based functions.)
+--
+--
+-- Test getting the distance between two points using earth_distance.
+--
+select earth_distance(ll_to_earth(0,0),ll_to_earth(0,0))::numeric(20,5);
+ earth_distance 
+----------------
+        0.00000
+(1 row)
+
+select earth_distance(ll_to_earth(0,0),ll_to_earth(0,180))::numeric(20,5);
+ earth_distance 
+----------------
+ 20037605.73216
+(1 row)
+
+select earth_distance(ll_to_earth(0,0),ll_to_earth(90,0))::numeric(20,5);
+ earth_distance 
+----------------
+ 10018802.86608
+(1 row)
+
+select earth_distance(ll_to_earth(0,0),ll_to_earth(0,90))::numeric(20,5);
+ earth_distance 
+----------------
+ 10018802.86608
+(1 row)
+
+select earth_distance(ll_to_earth(0,0),ll_to_earth(0,1))::numeric(20,5);
+ earth_distance 
+----------------
+   111320.03185
+(1 row)
+
+select earth_distance(ll_to_earth(0,0),ll_to_earth(1,0))::numeric(20,5);
+ earth_distance 
+----------------
+   111320.03185
+(1 row)
+
+select earth_distance(ll_to_earth(30,0),ll_to_earth(30,1))::numeric(20,5);
+ earth_distance 
+----------------
+    96405.66962
+(1 row)
+
+select earth_distance(ll_to_earth(30,0),ll_to_earth(31,0))::numeric(20,5);
+ earth_distance 
+----------------
+   111320.03185
+(1 row)
+
+select earth_distance(ll_to_earth(60,0),ll_to_earth(60,1))::numeric(20,5);
+ earth_distance 
+----------------
+    55659.48608
+(1 row)
+
+select earth_distance(ll_to_earth(60,0),ll_to_earth(61,0))::numeric(20,5);
+ earth_distance 
+----------------
+   111320.03185
+(1 row)
+
+select earth_distance(ll_to_earth(41.8,87.6),ll_to_earth(35.1,106.7))::numeric(20,5);
+ earth_distance 
+----------------
+  1819303.21265
+(1 row)
+
+select (earth_distance(ll_to_earth(41.8,87.6),ll_to_earth(35.1,106.7))*
+      100./2.54/12./5280.)::numeric(20,5);
+  numeric   
+------------
+ 1130.46261
+(1 row)
+
 --
 -- Test getting the distance between two points using geo_distance.
 --
-SELECT geo_distance('(0,0)'::point,'(0,0)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(0,0)'::point)::numeric(20,5);
  geo_distance 
 --------------
       0.00000
 (1 row)
 
-SELECT geo_distance('(0,0)'::point,'(180,0)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(180,0)'::point)::numeric(20,5);
  geo_distance 
 --------------
   12436.77274
 (1 row)
 
-SELECT geo_distance('(0,0)'::point,'(0,90)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(0,90)'::point)::numeric(20,5);
  geo_distance 
 --------------
    6218.38637
 (1 row)
 
-SELECT geo_distance('(0,0)'::point,'(90,0)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(90,0)'::point)::numeric(20,5);
  geo_distance 
 --------------
    6218.38637
 (1 row)
 
-SELECT geo_distance('(0,0)'::point,'(1,0)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(1,0)'::point)::numeric(20,5);
  geo_distance 
 --------------
      69.09318
 (1 row)
 
-SELECT geo_distance('(0,0)'::point,'(0,1)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(0,1)'::point)::numeric(20,5);
  geo_distance 
 --------------
      69.09318
 (1 row)
 
-SELECT geo_distance('(0,30)'::point,'(1,30)'::point)::numeric(20,5);
+select geo_distance('(0,30)'::point,'(1,30)'::point)::numeric(20,5);
  geo_distance 
 --------------
      59.83626
 (1 row)
 
-SELECT geo_distance('(0,30)'::point,'(0,31)'::point)::numeric(20,5);
+select geo_distance('(0,30)'::point,'(0,31)'::point)::numeric(20,5);
  geo_distance 
 --------------
      69.09318
 (1 row)
 
-SELECT geo_distance('(0,60)'::point,'(1,60)'::point)::numeric(20,5);
+select geo_distance('(0,60)'::point,'(1,60)'::point)::numeric(20,5);
  geo_distance 
 --------------
      34.54626
 (1 row)
 
-SELECT geo_distance('(0,60)'::point,'(0,61)'::point)::numeric(20,5);
+select geo_distance('(0,60)'::point,'(0,61)'::point)::numeric(20,5);
  geo_distance 
 --------------
      69.09318
 (1 row)
 
-SELECT geo_distance('(87.6,41.8)'::point,'(106.7,35.1)'::point)::numeric(20,5);
+select geo_distance('(87.6,41.8)'::point,'(106.7,35.1)'::point)::numeric(20,5);
  geo_distance 
 --------------
    1129.18983
 (1 row)
 
-SELECT (geo_distance('(87.6,41.8)'::point,'(106.7,35.1)'::point)*5280.*12.*2.54/100.)::numeric(20,5);
+select (geo_distance('(87.6,41.8)'::point,'(106.7,35.1)'::point)*5280.*12.*2.54/100.)::numeric(20,5);
     numeric    
 ---------------
  1817254.87730
@@ -84,75 +631,318 @@ SELECT (geo_distance('(87.6,41.8)'::point,'(106.7,35.1)'::point)*5280.*12.*2.54/
 --
 -- Test getting the distance between two points using the <@> operator.
 --
-SELECT ('(0,0)'::point <@> '(0,0)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(0,0)'::point)::numeric(20,5);
  numeric 
 ---------
  0.00000
 (1 row)
 
-SELECT ('(0,0)'::point <@> '(180,0)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(180,0)'::point)::numeric(20,5);
    numeric   
 -------------
  12436.77274
 (1 row)
 
-SELECT ('(0,0)'::point <@> '(0,90)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(0,90)'::point)::numeric(20,5);
   numeric   
 ------------
  6218.38637
 (1 row)
 
-SELECT ('(0,0)'::point <@> '(90,0)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(90,0)'::point)::numeric(20,5);
   numeric   
 ------------
  6218.38637
 (1 row)
 
-SELECT ('(0,0)'::point <@> '(1,0)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(1,0)'::point)::numeric(20,5);
  numeric  
 ----------
  69.09318
 (1 row)
 
-SELECT ('(0,0)'::point <@> '(0,1)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(0,1)'::point)::numeric(20,5);
  numeric  
 ----------
  69.09318
 (1 row)
 
-SELECT ('(0,30)'::point <@> '(1,30)'::point)::numeric(20,5);
+select ('(0,30)'::point <@> '(1,30)'::point)::numeric(20,5);
  numeric  
 ----------
  59.83626
 (1 row)
 
-SELECT ('(0,30)'::point <@> '(0,31)'::point)::numeric(20,5);
+select ('(0,30)'::point <@> '(0,31)'::point)::numeric(20,5);
  numeric  
 ----------
  69.09318
 (1 row)
 
-SELECT ('(0,60)'::point <@> '(1,60)'::point)::numeric(20,5);
+select ('(0,60)'::point <@> '(1,60)'::point)::numeric(20,5);
  numeric  
 ----------
  34.54626
 (1 row)
 
-SELECT ('(0,60)'::point <@> '(0,61)'::point)::numeric(20,5);
+select ('(0,60)'::point <@> '(0,61)'::point)::numeric(20,5);
  numeric  
 ----------
  69.09318
 (1 row)
 
-SELECT ('(87.6,41.8)'::point <@> '(106.7,35.1)'::point)::numeric(20,5);
+select ('(87.6,41.8)'::point <@> '(106.7,35.1)'::point)::numeric(20,5);
   numeric   
 ------------
  1129.18983
 (1 row)
 
-SELECT (('(87.6,41.8)'::point <@> '(106.7,35.1)'::point)*5280.*12.*2.54/100.)::numeric(20,5);
+select (('(87.6,41.8)'::point <@> '(106.7,35.1)'::point)*5280.*12.*2.54/100.)::numeric(20,5);
     numeric    
 ---------------
  1817254.87730
 (1 row)
 
+--
+-- Test getting a bounding box around points.
+--
+select cube_ll_coord(earth_box(ll_to_earth(0,0),112000),1)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),112000),2)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),112000),3)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),112000),1)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),112000),2)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),112000),3)::numeric(20,5);
+ cube_ll_coord | cube_ll_coord | cube_ll_coord | cube_ur_coord | cube_ur_coord | cube_ur_coord 
+---------------+---------------+---------------+---------------+---------------+---------------
+ 6266169.43896 | -111998.56104 | -111998.56104 | 6490166.56104 |  111998.56104 |  111998.56104
+(1 row)
+
+select cube_ll_coord(earth_box(ll_to_earth(0,0),pi()*earth()),1)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),pi()*earth()),2)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),pi()*earth()),3)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),pi()*earth()),1)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),pi()*earth()),2)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),pi()*earth()),3)::numeric(20,5);
+ cube_ll_coord  |  cube_ll_coord  |  cube_ll_coord  | cube_ur_coord  | cube_ur_coord  | cube_ur_coord  
+----------------+-----------------+-----------------+----------------+----------------+----------------
+ -6378168.00000 | -12756336.00000 | -12756336.00000 | 19134504.00000 | 12756336.00000 | 12756336.00000
+(1 row)
+
+select cube_ll_coord(earth_box(ll_to_earth(0,0),10*earth()),1)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),10*earth()),2)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),10*earth()),3)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),10*earth()),1)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),10*earth()),2)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),10*earth()),3)::numeric(20,5);
+ cube_ll_coord  |  cube_ll_coord  |  cube_ll_coord  | cube_ur_coord  | cube_ur_coord  | cube_ur_coord  
+----------------+-----------------+-----------------+----------------+----------------+----------------
+ -6378168.00000 | -12756336.00000 | -12756336.00000 | 19134504.00000 | 12756336.00000 | 12756336.00000
+(1 row)
+
+--
+-- Test for points that should be in bounding boxes.
+--
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,1))*1.00001) @
+       ll_to_earth(0,1);
+ ?column? 
+----------
+ t
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.1))*1.00001) @
+       ll_to_earth(0,0.1);
+ ?column? 
+----------
+ t
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.01))*1.00001) @
+       ll_to_earth(0,0.01);
+ ?column? 
+----------
+ t
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.001))*1.00001) @
+       ll_to_earth(0,0.001);
+ ?column? 
+----------
+ t
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.0001))*1.00001) @
+       ll_to_earth(0,0.0001);
+ ?column? 
+----------
+ t
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0.0001,0.0001))*1.00001) @
+       ll_to_earth(0.0001,0.0001);
+ ?column? 
+----------
+ t
+(1 row)
+
+select earth_box(ll_to_earth(45,45),
+       earth_distance(ll_to_earth(45,45),ll_to_earth(45.0001,45.0001))*1.00001) @
+       ll_to_earth(45.0001,45.0001);
+ ?column? 
+----------
+ t
+(1 row)
+
+select earth_box(ll_to_earth(90,180),
+       earth_distance(ll_to_earth(90,180),ll_to_earth(90.0001,180.0001))*1.00001) @
+       ll_to_earth(90.0001,180.0001);
+ ?column? 
+----------
+ t
+(1 row)
+
+--
+-- Test for points that shouldn't be in bounding boxes. Note that we need
+-- to make points way outside, since some points close may be in the box
+-- but further away than the distance we are testing.
+--
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,1))*.57735) @
+       ll_to_earth(0,1);
+ ?column? 
+----------
+ f
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.1))*.57735) @
+       ll_to_earth(0,0.1);
+ ?column? 
+----------
+ f
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.01))*.57735) @
+       ll_to_earth(0,0.01);
+ ?column? 
+----------
+ f
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.001))*.57735) @
+       ll_to_earth(0,0.001);
+ ?column? 
+----------
+ f
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.0001))*.57735) @
+       ll_to_earth(0,0.0001);
+ ?column? 
+----------
+ f
+(1 row)
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0.0001,0.0001))*.57735) @
+       ll_to_earth(0.0001,0.0001);
+ ?column? 
+----------
+ f
+(1 row)
+
+select earth_box(ll_to_earth(45,45),
+       earth_distance(ll_to_earth(45,45),ll_to_earth(45.0001,45.0001))*.57735) @
+       ll_to_earth(45.0001,45.0001);
+ ?column? 
+----------
+ f
+(1 row)
+
+select earth_box(ll_to_earth(90,180),
+       earth_distance(ll_to_earth(90,180),ll_to_earth(90.0001,180.0001))*.57735) @
+       ll_to_earth(90.0001,180.0001);
+ ?column? 
+----------
+ f
+(1 row)
+
+--
+-- Test the recommended constraints.
+--
+select is_point(ll_to_earth(0,0));
+ERROR:  Function is_point(earth) does not exist
+       Unable to identify a function that satisfies the given argument types
+       You may need to add explicit typecasts
+select cube_dim(ll_to_earth(0,0)) <= 3;
+ ?column? 
+----------
+ t
+(1 row)
+
+select abs(cube_distance(ll_to_earth(0,0), '(0)'::cube) / earth() - 1) <
+       '10e-12'::float8;
+ ?column? 
+----------
+ t
+(1 row)
+
+select is_point(ll_to_earth(30,60));
+ERROR:  Function is_point(earth) does not exist
+       Unable to identify a function that satisfies the given argument types
+       You may need to add explicit typecasts
+select cube_dim(ll_to_earth(30,60)) <= 3;
+ ?column? 
+----------
+ t
+(1 row)
+
+select abs(cube_distance(ll_to_earth(30,60), '(0)'::cube) / earth() - 1) <
+       '10e-12'::float8;
+ ?column? 
+----------
+ t
+(1 row)
+
+select is_point(ll_to_earth(60,90));
+ERROR:  Function is_point(earth) does not exist
+       Unable to identify a function that satisfies the given argument types
+       You may need to add explicit typecasts
+select cube_dim(ll_to_earth(60,90)) <= 3;
+ ?column? 
+----------
+ t
+(1 row)
+
+select abs(cube_distance(ll_to_earth(60,90), '(0)'::cube) / earth() - 1) <
+       '10e-12'::float8;
+ ?column? 
+----------
+ t
+(1 row)
+
+select is_point(ll_to_earth(-30,-90));
+ERROR:  Function is_point(earth) does not exist
+       Unable to identify a function that satisfies the given argument types
+       You may need to add explicit typecasts
+select cube_dim(ll_to_earth(-30,-90)) <= 3;
+ ?column? 
+----------
+ t
+(1 row)
+
+select abs(cube_distance(ll_to_earth(-30,-90), '(0)'::cube) / earth() - 1) <
+       '10e-12'::float8;
+ ?column? 
+----------
+ t
+(1 row)
+
index b3b88f701dd72bbdfde4cfa61d1f6858946d7eab..29f3c47136502f97f3d7f20824596f1dacc4fa73 100644 (file)
 
 --
 -- first, define the datatype.  Turn off echoing so that expected file
--- does not depend on contents of earthdistance.sql.
+-- does not depend on contents of earthdistance.sql or cube.sql.
 --
 \set ECHO none
+\i ../cube/cube.sql
 \i earthdistance.sql
 \set ECHO all
 
+--
+-- The radius of the Earth we are using.
+--
+
+select earth()::numeric(20,5);
+
+--
+-- Convert straight line distances to great circle distances.
+--
+select (pi()*earth())::numeric(20,5);
+select sec_to_gc(0)::numeric(20,5);
+select sec_to_gc(2*earth())::numeric(20,5);
+select sec_to_gc(10*earth())::numeric(20,5);
+select sec_to_gc(-earth())::numeric(20,5);
+select sec_to_gc(1000)::numeric(20,5);
+select sec_to_gc(10000)::numeric(20,5);
+select sec_to_gc(100000)::numeric(20,5);
+select sec_to_gc(1000000)::numeric(20,5);
+
+--
+-- Convert great circle distances to straight line distances.
+--
+
+select gc_to_sec(0)::numeric(20,5);
+select gc_to_sec(sec_to_gc(2*earth()))::numeric(20,5);
+select gc_to_sec(10*earth())::numeric(20,5);
+select gc_to_sec(pi()*earth())::numeric(20,5);
+select gc_to_sec(-1000)::numeric(20,5);
+select gc_to_sec(1000)::numeric(20,5);
+select gc_to_sec(10000)::numeric(20,5);
+select gc_to_sec(100000)::numeric(20,5);
+select gc_to_sec(1000000)::numeric(20,5);
+
+--
+-- Set coordinates using latitude and longitude.
+-- Extract each coordinate separately so we can round them.
+--
+
+select cube_ll_coord(ll_to_earth(0,0),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,0),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,0),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(360,360),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(360,360),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(360,360),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(180,180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(180,180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(180,180),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(180,360),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(180,360),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(180,360),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(-180,-360),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-180,-360),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-180,-360),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(0,180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,180),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(0,-180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,-180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(0,-180),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(90,0),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(90,0),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(90,0),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(90,180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(90,180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(90,180),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(-90,0),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-90,0),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-90,0),3)::numeric(20,5);
+select cube_ll_coord(ll_to_earth(-90,180),1)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-90,180),2)::numeric(20,5),
+ cube_ll_coord(ll_to_earth(-90,180),3)::numeric(20,5);
+
+--
+-- Test getting the latitude of a location.
+--
+
+select latitude(ll_to_earth(0,0))::numeric(20,10);
+select latitude(ll_to_earth(45,0))::numeric(20,10);
+select latitude(ll_to_earth(90,0))::numeric(20,10);
+select latitude(ll_to_earth(-45,0))::numeric(20,10);
+select latitude(ll_to_earth(-90,0))::numeric(20,10);
+select latitude(ll_to_earth(0,90))::numeric(20,10);
+select latitude(ll_to_earth(45,90))::numeric(20,10);
+select latitude(ll_to_earth(90,90))::numeric(20,10);
+select latitude(ll_to_earth(-45,90))::numeric(20,10);
+select latitude(ll_to_earth(-90,90))::numeric(20,10);
+select latitude(ll_to_earth(0,180))::numeric(20,10);
+select latitude(ll_to_earth(45,180))::numeric(20,10);
+select latitude(ll_to_earth(90,180))::numeric(20,10);
+select latitude(ll_to_earth(-45,180))::numeric(20,10);
+select latitude(ll_to_earth(-90,180))::numeric(20,10);
+select latitude(ll_to_earth(0,-90))::numeric(20,10);
+select latitude(ll_to_earth(45,-90))::numeric(20,10);
+select latitude(ll_to_earth(90,-90))::numeric(20,10);
+select latitude(ll_to_earth(-45,-90))::numeric(20,10);
+select latitude(ll_to_earth(-90,-90))::numeric(20,10);
+
+--
+-- Test getting the longitude of a location.
+--
+
+select longitude(ll_to_earth(0,0))::numeric(20,10);
+select longitude(ll_to_earth(45,0))::numeric(20,10);
+select longitude(ll_to_earth(90,0))::numeric(20,10);
+select longitude(ll_to_earth(-45,0))::numeric(20,10);
+select longitude(ll_to_earth(-90,0))::numeric(20,10);
+select longitude(ll_to_earth(0,90))::numeric(20,10);
+select longitude(ll_to_earth(45,90))::numeric(20,10);
+select longitude(ll_to_earth(90,90))::numeric(20,10);
+select longitude(ll_to_earth(-45,90))::numeric(20,10);
+select longitude(ll_to_earth(-90,90))::numeric(20,10);
+select longitude(ll_to_earth(0,180))::numeric(20,10);
+select longitude(ll_to_earth(45,180))::numeric(20,10);
+select longitude(ll_to_earth(90,180))::numeric(20,10);
+select longitude(ll_to_earth(-45,180))::numeric(20,10);
+select longitude(ll_to_earth(-90,180))::numeric(20,10);
+select longitude(ll_to_earth(0,-90))::numeric(20,10);
+select longitude(ll_to_earth(45,-90))::numeric(20,10);
+select longitude(ll_to_earth(90,-90))::numeric(20,10);
+select longitude(ll_to_earth(-45,-90))::numeric(20,10);
+select longitude(ll_to_earth(-90,-90))::numeric(20,10);
+
+--
+-- For the distance tests the following is some real life data.
+--
+-- Chicago has a latitude of 41.8 and a longitude of 87.6.
+-- Albuquerque has a latitude of 35.1 and a longitude of 106.7.
+-- (Note that latitude and longitude are specified differently
+-- in the cube based functions than for the point based functions.)
+--
+
+--
+-- Test getting the distance between two points using earth_distance.
+--
+
+select earth_distance(ll_to_earth(0,0),ll_to_earth(0,0))::numeric(20,5);
+select earth_distance(ll_to_earth(0,0),ll_to_earth(0,180))::numeric(20,5);
+select earth_distance(ll_to_earth(0,0),ll_to_earth(90,0))::numeric(20,5);
+select earth_distance(ll_to_earth(0,0),ll_to_earth(0,90))::numeric(20,5);
+select earth_distance(ll_to_earth(0,0),ll_to_earth(0,1))::numeric(20,5);
+select earth_distance(ll_to_earth(0,0),ll_to_earth(1,0))::numeric(20,5);
+select earth_distance(ll_to_earth(30,0),ll_to_earth(30,1))::numeric(20,5);
+select earth_distance(ll_to_earth(30,0),ll_to_earth(31,0))::numeric(20,5);
+select earth_distance(ll_to_earth(60,0),ll_to_earth(60,1))::numeric(20,5);
+select earth_distance(ll_to_earth(60,0),ll_to_earth(61,0))::numeric(20,5);
+select earth_distance(ll_to_earth(41.8,87.6),ll_to_earth(35.1,106.7))::numeric(20,5);
+select (earth_distance(ll_to_earth(41.8,87.6),ll_to_earth(35.1,106.7))*
+      100./2.54/12./5280.)::numeric(20,5);
+
 --
 -- Test getting the distance between two points using geo_distance.
 --
 
-SELECT geo_distance('(0,0)'::point,'(0,0)'::point)::numeric(20,5);
-SELECT geo_distance('(0,0)'::point,'(180,0)'::point)::numeric(20,5);
-SELECT geo_distance('(0,0)'::point,'(0,90)'::point)::numeric(20,5);
-SELECT geo_distance('(0,0)'::point,'(90,0)'::point)::numeric(20,5);
-SELECT geo_distance('(0,0)'::point,'(1,0)'::point)::numeric(20,5);
-SELECT geo_distance('(0,0)'::point,'(0,1)'::point)::numeric(20,5);
-SELECT geo_distance('(0,30)'::point,'(1,30)'::point)::numeric(20,5);
-SELECT geo_distance('(0,30)'::point,'(0,31)'::point)::numeric(20,5);
-SELECT geo_distance('(0,60)'::point,'(1,60)'::point)::numeric(20,5);
-SELECT geo_distance('(0,60)'::point,'(0,61)'::point)::numeric(20,5);
-SELECT geo_distance('(87.6,41.8)'::point,'(106.7,35.1)'::point)::numeric(20,5);
-SELECT (geo_distance('(87.6,41.8)'::point,'(106.7,35.1)'::point)*5280.*12.*2.54/100.)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(0,0)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(180,0)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(0,90)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(90,0)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(1,0)'::point)::numeric(20,5);
+select geo_distance('(0,0)'::point,'(0,1)'::point)::numeric(20,5);
+select geo_distance('(0,30)'::point,'(1,30)'::point)::numeric(20,5);
+select geo_distance('(0,30)'::point,'(0,31)'::point)::numeric(20,5);
+select geo_distance('(0,60)'::point,'(1,60)'::point)::numeric(20,5);
+select geo_distance('(0,60)'::point,'(0,61)'::point)::numeric(20,5);
+select geo_distance('(87.6,41.8)'::point,'(106.7,35.1)'::point)::numeric(20,5);
+select (geo_distance('(87.6,41.8)'::point,'(106.7,35.1)'::point)*5280.*12.*2.54/100.)::numeric(20,5);
 
 --
 -- Test getting the distance between two points using the <@> operator.
 --
 
-SELECT ('(0,0)'::point <@> '(0,0)'::point)::numeric(20,5);
-SELECT ('(0,0)'::point <@> '(180,0)'::point)::numeric(20,5);
-SELECT ('(0,0)'::point <@> '(0,90)'::point)::numeric(20,5);
-SELECT ('(0,0)'::point <@> '(90,0)'::point)::numeric(20,5);
-SELECT ('(0,0)'::point <@> '(1,0)'::point)::numeric(20,5);
-SELECT ('(0,0)'::point <@> '(0,1)'::point)::numeric(20,5);
-SELECT ('(0,30)'::point <@> '(1,30)'::point)::numeric(20,5);
-SELECT ('(0,30)'::point <@> '(0,31)'::point)::numeric(20,5);
-SELECT ('(0,60)'::point <@> '(1,60)'::point)::numeric(20,5);
-SELECT ('(0,60)'::point <@> '(0,61)'::point)::numeric(20,5);
-SELECT ('(87.6,41.8)'::point <@> '(106.7,35.1)'::point)::numeric(20,5);
-SELECT (('(87.6,41.8)'::point <@> '(106.7,35.1)'::point)*5280.*12.*2.54/100.)::numeric(20,5);
+select ('(0,0)'::point <@> '(0,0)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(180,0)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(0,90)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(90,0)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(1,0)'::point)::numeric(20,5);
+select ('(0,0)'::point <@> '(0,1)'::point)::numeric(20,5);
+select ('(0,30)'::point <@> '(1,30)'::point)::numeric(20,5);
+select ('(0,30)'::point <@> '(0,31)'::point)::numeric(20,5);
+select ('(0,60)'::point <@> '(1,60)'::point)::numeric(20,5);
+select ('(0,60)'::point <@> '(0,61)'::point)::numeric(20,5);
+select ('(87.6,41.8)'::point <@> '(106.7,35.1)'::point)::numeric(20,5);
+select (('(87.6,41.8)'::point <@> '(106.7,35.1)'::point)*5280.*12.*2.54/100.)::numeric(20,5);
+
+--
+-- Test getting a bounding box around points.
+--
+
+select cube_ll_coord(earth_box(ll_to_earth(0,0),112000),1)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),112000),2)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),112000),3)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),112000),1)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),112000),2)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),112000),3)::numeric(20,5);
+select cube_ll_coord(earth_box(ll_to_earth(0,0),pi()*earth()),1)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),pi()*earth()),2)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),pi()*earth()),3)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),pi()*earth()),1)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),pi()*earth()),2)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),pi()*earth()),3)::numeric(20,5);
+select cube_ll_coord(earth_box(ll_to_earth(0,0),10*earth()),1)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),10*earth()),2)::numeric(20,5),
+       cube_ll_coord(earth_box(ll_to_earth(0,0),10*earth()),3)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),10*earth()),1)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),10*earth()),2)::numeric(20,5),
+       cube_ur_coord(earth_box(ll_to_earth(0,0),10*earth()),3)::numeric(20,5);
+
+--
+-- Test for points that should be in bounding boxes.
+--
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,1))*1.00001) @
+       ll_to_earth(0,1);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.1))*1.00001) @
+       ll_to_earth(0,0.1);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.01))*1.00001) @
+       ll_to_earth(0,0.01);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.001))*1.00001) @
+       ll_to_earth(0,0.001);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.0001))*1.00001) @
+       ll_to_earth(0,0.0001);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0.0001,0.0001))*1.00001) @
+       ll_to_earth(0.0001,0.0001);
+select earth_box(ll_to_earth(45,45),
+       earth_distance(ll_to_earth(45,45),ll_to_earth(45.0001,45.0001))*1.00001) @
+       ll_to_earth(45.0001,45.0001);
+select earth_box(ll_to_earth(90,180),
+       earth_distance(ll_to_earth(90,180),ll_to_earth(90.0001,180.0001))*1.00001) @
+       ll_to_earth(90.0001,180.0001);
+
+--
+-- Test for points that shouldn't be in bounding boxes. Note that we need
+-- to make points way outside, since some points close may be in the box
+-- but further away than the distance we are testing.
+--
+
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,1))*.57735) @
+       ll_to_earth(0,1);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.1))*.57735) @
+       ll_to_earth(0,0.1);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.01))*.57735) @
+       ll_to_earth(0,0.01);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.001))*.57735) @
+       ll_to_earth(0,0.001);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0,0.0001))*.57735) @
+       ll_to_earth(0,0.0001);
+select earth_box(ll_to_earth(0,0),
+       earth_distance(ll_to_earth(0,0),ll_to_earth(0.0001,0.0001))*.57735) @
+       ll_to_earth(0.0001,0.0001);
+select earth_box(ll_to_earth(45,45),
+       earth_distance(ll_to_earth(45,45),ll_to_earth(45.0001,45.0001))*.57735) @
+       ll_to_earth(45.0001,45.0001);
+select earth_box(ll_to_earth(90,180),
+       earth_distance(ll_to_earth(90,180),ll_to_earth(90.0001,180.0001))*.57735) @
+       ll_to_earth(90.0001,180.0001);
+
+--
+-- Test the recommended constraints.
+--
+
+select is_point(ll_to_earth(0,0));
+select cube_dim(ll_to_earth(0,0)) <= 3;
+select abs(cube_distance(ll_to_earth(0,0), '(0)'::cube) / earth() - 1) <
+       '10e-12'::float8;
+select is_point(ll_to_earth(30,60));
+select cube_dim(ll_to_earth(30,60)) <= 3;
+select abs(cube_distance(ll_to_earth(30,60), '(0)'::cube) / earth() - 1) <
+       '10e-12'::float8;
+select is_point(ll_to_earth(60,90));
+select cube_dim(ll_to_earth(60,90)) <= 3;
+select abs(cube_distance(ll_to_earth(60,90), '(0)'::cube) / earth() - 1) <
+       '10e-12'::float8;
+select is_point(ll_to_earth(-30,-90));
+select cube_dim(ll_to_earth(-30,-90)) <= 3;
+select abs(cube_distance(ll_to_earth(-30,-90), '(0)'::cube) / earth() - 1) <
+       '10e-12'::float8;