PostGIS  2.2.8dev-r@@SVN_REVISION@@
lwgeom_spheroid.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * Copyright (C) 2001-2003 Refractions Research Inc.
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include "postgres.h"
14 
15 
16 #include <math.h>
17 #include <float.h>
18 #include <string.h>
19 #include <stdio.h>
20 #include <errno.h>
21 
22 #include "access/gist.h"
23 #include "access/itup.h"
24 
25 #include "fmgr.h"
26 #include "utils/elog.h"
27 
28 #include "../postgis_config.h"
29 #include "liblwgeom.h"
30 #include "lwgeom_pg.h"
31 
32 #define SHOW_DIGS_DOUBLE 15
33 #define MAX_DIGS_DOUBLE (SHOW_DIGS_DOUBLE + 6 + 1 + 3 +1)
34 
35 /*
36  * distance from -126 49 to -126 49.011096139863
37  * in 'SPHEROID["GRS_1980",6378137,298.257222101]'
38  * is 1234.000
39  */
40 
41 
42 /* PG-exposed */
43 Datum ellipsoid_in(PG_FUNCTION_ARGS);
44 Datum ellipsoid_out(PG_FUNCTION_ARGS);
45 Datum LWGEOM_length2d_ellipsoid(PG_FUNCTION_ARGS);
46 Datum LWGEOM_length_ellipsoid_linestring(PG_FUNCTION_ARGS);
47 Datum LWGEOM_distance_ellipsoid(PG_FUNCTION_ARGS);
48 Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS);
49 Datum geometry_distance_spheroid(PG_FUNCTION_ARGS);
50 
51 /* internal */
52 double distance_sphere_method(double lat1, double long1,double lat2,double long2, SPHEROID *sphere);
53 double distance_ellipse_calculation(double lat1, double long1, double lat2, double long2, SPHEROID *sphere);
54 double distance_ellipse(double lat1, double long1, double lat2, double long2, SPHEROID *sphere);
55 double deltaLongitude(double azimuth, double sigma, double tsm,SPHEROID *sphere);
56 double mu2(double azimuth,SPHEROID *sphere);
57 double bigA(double u2);
58 double bigB(double u2);
59 
60 
61 /*
62  * Use the WKT definition of an ellipsoid
63  * ie. SPHEROID["name",A,rf] or SPHEROID("name",A,rf)
64  * SPHEROID["GRS_1980",6378137,298.257222101]
65  * wkt says you can use "(" or "["
66  */
68 Datum ellipsoid_in(PG_FUNCTION_ARGS)
69 {
70  char *str = PG_GETARG_CSTRING(0);
71  SPHEROID *sphere = (SPHEROID *) palloc(sizeof(SPHEROID));
72  int nitems;
73  double rf;
74 
75  memset(sphere,0, sizeof(SPHEROID));
76 
77  if (strstr(str,"SPHEROID") != str )
78  {
79  elog(ERROR,"SPHEROID parser - doesn't start with SPHEROID");
80  pfree(sphere);
81  PG_RETURN_NULL();
82  }
83 
84  nitems = sscanf(str,"SPHEROID[\"%19[^\"]\",%lf,%lf]",
85  sphere->name, &sphere->a, &rf);
86 
87  if ( nitems==0)
88  nitems = sscanf(str,"SPHEROID(\"%19[^\"]\",%lf,%lf)",
89  sphere->name, &sphere->a, &rf);
90 
91  if (nitems != 3)
92  {
93  elog(ERROR,"SPHEROID parser - couldnt parse the spheroid");
94  pfree(sphere);
95  PG_RETURN_NULL();
96  }
97 
98  sphere->f = 1.0/rf;
99  sphere->b = sphere->a - (1.0/rf)*sphere->a;
100  sphere->e_sq = ((sphere->a*sphere->a) - (sphere->b*sphere->b)) /
101  (sphere->a*sphere->a);
102  sphere->e = sqrt(sphere->e_sq);
103 
104  PG_RETURN_POINTER(sphere);
105 
106 }
107 
109 Datum ellipsoid_out(PG_FUNCTION_ARGS)
110 {
111  SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(0);
112  char *result;
113 
114  result = palloc(MAX_DIGS_DOUBLE + MAX_DIGS_DOUBLE + 20 + 9 + 2);
115 
116  sprintf(result,"SPHEROID(\"%s\",%.15g,%.15g)",
117  sphere->name, sphere->a, 1.0/sphere->f);
118 
119  PG_RETURN_CSTRING(result);
120 }
121 
122 /*
123  * support function for distance calc
124  * code is taken from David Skea
125  * Geographic Data BC, Province of British Columbia, Canada.
126  * Thanks to GDBC and David Skea for allowing this to be
127  * put in PostGIS.
128  */
129 double
130 deltaLongitude(double azimuth, double sigma, double tsm,SPHEROID *sphere)
131 {
132  /* compute the expansion C */
133  double das,C;
134  double ctsm,DL;
135 
136  das = cos(azimuth)*cos(azimuth);
137  C = sphere->f/16.0 * das * (4.0 + sphere->f * (4.0 - 3.0 * das));
138 
139  /* compute the difference in longitude */
140  ctsm = cos(tsm);
141  DL = ctsm + C * cos(sigma) * (-1.0 + 2.0 * ctsm*ctsm);
142  DL = sigma + C * sin(sigma) * DL;
143  return (1.0 - C) * sphere->f * sin(azimuth) * DL;
144 }
145 
146 
147 /*
148  * support function for distance calc
149  * code is taken from David Skea
150  * Geographic Data BC, Province of British Columbia, Canada.
151  * Thanks to GDBC and David Skea for allowing this to be
152  * put in PostGIS.
153  */
154 double
155 mu2(double azimuth,SPHEROID *sphere)
156 {
157  double e2;
158 
159  e2 = sqrt(sphere->a*sphere->a-sphere->b*sphere->b)/sphere->b;
160  return cos(azimuth)*cos(azimuth) * e2*e2;
161 }
162 
163 
164 /*
165  * Support function for distance calc
166  * code is taken from David Skea
167  * Geographic Data BC, Province of British Columbia, Canada.
168  * Thanks to GDBC and David Skea for allowing this to be
169  * put in PostGIS.
170  */
171 double
172 bigA(double u2)
173 {
174  return 1.0 + u2/256.0 * (64.0 + u2 * (-12.0 + 5.0 * u2));
175 }
176 
177 
178 /*
179  * Support function for distance calc
180  * code is taken from David Skea
181  * Geographic Data BC, Province of British Columbia, Canada.
182  * Thanks to GDBC and David Skea for allowing this to be
183  * put in PostGIS.
184  */
185 double
186 bigB(double u2)
187 {
188  return u2/512.0 * (128.0 + u2 * (-64.0 + 37.0 * u2));
189 }
190 
191 
192 
193 double
194 distance_ellipse(double lat1, double long1,
195  double lat2, double long2, SPHEROID *sphere)
196 {
197  double result = 0;
198 #if POSTGIS_DEBUG_LEVEL >= 4
199  double result2 = 0;
200 #endif
201 
202  if ( (lat1==lat2) && (long1 == long2) )
203  {
204  return 0.0; /* same point, therefore zero distance */
205  }
206 
207  result = distance_ellipse_calculation(lat1,long1,lat2,long2,sphere);
208 
209 #if POSTGIS_DEBUG_LEVEL >= 4
210  result2 = distance_sphere_method(lat1, long1,lat2,long2, sphere);
211 
212  POSTGIS_DEBUGF(4, "delta = %lf, skae says: %.15lf,2 circle says: %.15lf",
213  (result2-result),result,result2);
214  POSTGIS_DEBUGF(4, "2 circle says: %.15lf",result2);
215 #endif
216 
217  if (result != result) /* NaN check
218  * (x==x for all x except NaN by IEEE definition)
219  */
220  {
221  result = distance_sphere_method(lat1, long1,
222  lat2,long2, sphere);
223  }
224 
225  return result;
226 }
227 
228 /*
229  * Given 2 lat/longs and ellipse, find the distance
230  * note original r = 1st, s=2nd location
231  */
232 double
233 distance_ellipse_calculation(double lat1, double long1,
234  double lat2, double long2, SPHEROID *sphere)
235 {
236  /*
237  * Code is taken from David Skea
238  * Geographic Data BC, Province of British Columbia, Canada.
239  * Thanks to GDBC and David Skea for allowing this to be
240  * put in PostGIS.
241  */
242 
243  double L1,L2,sinU1,sinU2,cosU1,cosU2;
244  double dl,dl1,dl2,dl3,cosdl1,sindl1;
245  double cosSigma,sigma,azimuthEQ,tsm;
246  double u2,A,B;
247  double dsigma;
248 
249  double TEMP;
250 
251  int iterations;
252 
253 
254  L1 = atan((1.0 - sphere->f ) * tan( lat1) );
255  L2 = atan((1.0 - sphere->f ) * tan( lat2) );
256  sinU1 = sin(L1);
257  sinU2 = sin(L2);
258  cosU1 = cos(L1);
259  cosU2 = cos(L2);
260 
261  dl = long2- long1;
262  dl1 = dl;
263  cosdl1 = cos(dl);
264  sindl1 = sin(dl);
265  iterations = 0;
266  do
267  {
268  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1;
269  sigma = acos(cosSigma);
270  azimuthEQ = asin((cosU1 * cosU2 * sindl1)/sin(sigma));
271 
272  /*
273  * Patch from Patrica Tozer to handle minor
274  * mathematical stability problem
275  */
276  TEMP = cosSigma - (2.0 * sinU1 * sinU2)/
277  (cos(azimuthEQ)*cos(azimuthEQ));
278  if (TEMP > 1)
279  {
280  TEMP = 1;
281  }
282  else if (TEMP < -1)
283  {
284  TEMP = -1;
285  }
286  tsm = acos(TEMP);
287 
288 
289  /* (old code?)
290  tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ)));
291  */
292 
293  dl2 = deltaLongitude(azimuthEQ, sigma, tsm,sphere);
294  dl3 = dl1 - (dl + dl2);
295  dl1 = dl + dl2;
296  cosdl1 = cos(dl1);
297  sindl1 = sin(dl1);
298  iterations++;
299  }
300  while ( (iterations<999) && (fabs(dl3) > 1.0e-032));
301 
302  /* compute expansions A and B */
303  u2 = mu2(azimuthEQ,sphere);
304  A = bigA(u2);
305  B = bigB(u2);
306 
307  /* compute length of geodesic */
308  dsigma = B * sin(sigma) * (cos(tsm) +
309  (B*cosSigma*(-1.0 + 2.0 * (cos(tsm)*cos(tsm))))/4.0);
310  return sphere->b * (A * (sigma - dsigma));
311 }
312 
313 
314 /*
315  * Find the "length of a geometry"
316  * length2d_spheroid(point, sphere) = 0
317  * length2d_spheroid(line, sphere) = length of line
318  * length2d_spheroid(polygon, sphere) = 0
319  * -- could make sense to return sum(ring perimeter)
320  * uses ellipsoidal math to find the distance
321  * x's are longitude, and y's are latitude - both in decimal degrees
322  */
324 Datum LWGEOM_length2d_ellipsoid(PG_FUNCTION_ARGS)
325 {
326  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
327  SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1);
328  LWGEOM *lwgeom = lwgeom_from_gserialized(geom);
329  double dist = lwgeom_length_spheroid(lwgeom, sphere);
330  lwgeom_free(lwgeom);
331  PG_FREE_IF_COPY(geom, 0);
332  PG_RETURN_FLOAT8(dist);
333 }
334 
335 
336 /*
337  * Find the "length of a geometry"
338  *
339  * length2d_spheroid(point, sphere) = 0
340  * length2d_spheroid(line, sphere) = length of line
341  * length2d_spheroid(polygon, sphere) = 0
342  * -- could make sense to return sum(ring perimeter)
343  * uses ellipsoidal math to find the distance
344  * x's are longitude, and y's are latitude - both in decimal degrees
345  */
347 Datum LWGEOM_length_ellipsoid_linestring(PG_FUNCTION_ARGS)
348 {
349  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
350  LWGEOM *lwgeom = lwgeom_from_gserialized(geom);
351  SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1);
352  double length = 0.0;
353 
354  /* EMPTY things have no length */
355  if ( lwgeom_is_empty(lwgeom) )
356  {
357  lwgeom_free(lwgeom);
358  PG_RETURN_FLOAT8(0.0);
359  }
360 
361  length = lwgeom_length_spheroid(lwgeom, sphere);
362  lwgeom_free(lwgeom);
363  PG_FREE_IF_COPY(geom, 0);
364 
365  /* Something went wrong... */
366  if ( length < 0.0 )
367  {
368  elog(ERROR, "lwgeom_length_spheroid returned length < 0.0");
369  PG_RETURN_NULL();
370  }
371 
372  /* Clean up */
373  PG_RETURN_FLOAT8(length);
374 }
375 
376 /*
377  * For some lat/long points, the above method doesn't calculate the distance very well.
378  * Typically this is for two lat/long points that are very very close together (<10cm).
379  * This gets worse closer to the equator.
380  *
381  * This method works very well for very close together points, not so well if they're
382  * far away (>1km).
383  *
384  * METHOD:
385  * We create two circles (with Radius R and Radius S) and use these to calculate the distance.
386  *
387  * The first (R) is basically a (north-south) line of longitude.
388  * Its radius is approximated by looking at the ellipse. Near the equator R = 'a' (earth's major axis)
389  * near the pole R = 'b' (earth's minor axis).
390  *
391  * The second (S) is basically a (east-west) line of lattitude.
392  * Its radius runs from 'a' (major axis) at the equator, and near 0 at the poles.
393  *
394  *
395  * North pole
396  * *
397  * *
398  * *\--S--
399  * * R +
400  * * \ +
401  * * A\ +
402  * * ------ \ Equator/centre of earth
403  * *
404  * *
405  * *
406  * *
407  * *
408  * *
409  * South pole
410  * (side view of earth)
411  *
412  * Angle A is lat1
413  * R is the distance from the centre of the earth to the lat1/long1 point on the surface
414  * of the Earth.
415  * S is the circle-of-lattitude. Its calculated from the right triangle defined by
416  * the angle (90-A), and the hypothenus R.
417  *
418  *
419  *
420  * Once R and S have been calculated, the actual distance between the two points can be
421  * calculated.
422  *
423  * We dissolve the vector from lat1,long1 to lat2,long2 into its X and Y components (called DeltaX,DeltaY).
424  * The actual distance that these angle-based measurements represent is taken from the two
425  * circles we just calculated; R (for deltaY) and S (for deltaX).
426  *
427  * (if deltaX is 1 degrees, then that distance represents 1/360 of a circle of radius S.)
428  *
429  *
430  * Parts taken from PROJ4 - geodetic_to_geocentric() (for calculating Rn)
431  *
432  * remember that lat1/long1/lat2/long2 are comming in a *RADIANS* not degrees.
433  *
434  * By Patricia Tozer and Dave Blasby
435  *
436  * This is also called the "curvature method".
437  */
438 
439 double distance_sphere_method(double lat1, double long1,double lat2,double long2, SPHEROID *sphere)
440 {
441  double R,S,X,Y,deltaX,deltaY;
442 
443  double distance = 0.0;
444  double sin_lat = sin(lat1);
445  double sin2_lat = sin_lat * sin_lat;
446 
447  double Geocent_a = sphere->a;
448  double Geocent_e2 = sphere->e_sq;
449 
450  R = Geocent_a / (sqrt(1.0e0 - Geocent_e2 * sin2_lat));
451  /* 90 - lat1, but in radians */
452  S = R * sin(M_PI_2 - lat1) ;
453 
454  deltaX = long2 - long1; /* in rads */
455  deltaY = lat2 - lat1; /* in rads */
456 
457  /* think: a % of 2*pi*S */
458  X = deltaX/(2.0*M_PI) * 2 * M_PI * S;
459  Y = deltaY/(2.0*M_PI) * 2 * M_PI * R;
460 
461  distance = sqrt((X * X + Y * Y));
462 
463  return distance;
464 }
465 
466 /*
467  * distance (geometry,geometry, sphere)
468  */
470 Datum geometry_distance_spheroid(PG_FUNCTION_ARGS)
471 {
472  GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
473  GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
474  SPHEROID *sphere = (SPHEROID *)PG_GETARG_POINTER(2);
475  int type1 = gserialized_get_type(geom1);
476  int type2 = gserialized_get_type(geom2);
477  bool use_spheroid = PG_GETARG_BOOL(3);
478  LWGEOM *lwgeom1, *lwgeom2;
479  double distance;
480 
481  /* Calculate some other parameters on the spheroid */
482  spheroid_init(sphere, sphere->a, sphere->b);
483 
485 
486  /* Catch sphere special case and re-jig spheroid appropriately */
487  if ( ! use_spheroid )
488  {
489  sphere->a = sphere->b = sphere->radius;
490  }
491 
492  if ( ! ( type1 == POINTTYPE || type1 == LINETYPE || type1 == POLYGONTYPE ||
493  type1 == MULTIPOINTTYPE || type1 == MULTILINETYPE || type1 == MULTIPOLYGONTYPE ))
494  {
495  elog(ERROR, "geometry_distance_spheroid: Only point/line/polygon supported.\n");
496  PG_RETURN_NULL();
497  }
498 
499  if ( ! ( type2 == POINTTYPE || type2 == LINETYPE || type2 == POLYGONTYPE ||
500  type2 == MULTIPOINTTYPE || type2 == MULTILINETYPE || type2 == MULTIPOLYGONTYPE ))
501  {
502  elog(ERROR, "geometry_distance_spheroid: Only point/line/polygon supported.\n");
503  PG_RETURN_NULL();
504  }
505 
506  /* Get #LWGEOM structures */
507  lwgeom1 = lwgeom_from_gserialized(geom1);
508  lwgeom2 = lwgeom_from_gserialized(geom2);
509 
510  /* We are going to be calculating geodetic distances */
511  lwgeom_set_geodetic(lwgeom1, LW_TRUE);
512  lwgeom_set_geodetic(lwgeom2, LW_TRUE);
513 
514  distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, sphere, 0.0);
515 
516  PG_RETURN_FLOAT8(distance);
517 
518 }
519 
521 Datum LWGEOM_distance_ellipsoid(PG_FUNCTION_ARGS)
522 {
523  PG_RETURN_DATUM(DirectFunctionCall4(geometry_distance_spheroid,
524  PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), PG_GETARG_DATUM(2), BoolGetDatum(TRUE)));
525 }
526 
528 Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS)
529 {
530  SPHEROID s;
531 
532  /* Init to WGS84 */
533  spheroid_init(&s, 6378137.0, 6356752.314245179498);
534  s.a = s.b = s.radius;
535 
536  PG_RETURN_DATUM(DirectFunctionCall4(geometry_distance_spheroid,
537  PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), PointerGetDatum(&s), BoolGetDatum(FALSE)));
538 }
539 
#define LINETYPE
Definition: liblwgeom.h:71
uint32_t gserialized_get_type(const GSERIALIZED *s)
Extract the geometry type from the serialized form (it hides in the anonymous data area...
Definition: g_serialized.c:55
Datum LWGEOM_length_ellipsoid_linestring(PG_FUNCTION_ARGS)
void spheroid_init(SPHEROID *s, double a, double b)
Initialize a spheroid object for use in geodetic functions.
Definition: lwspheroid.c:25
#define MAX_DIGS_DOUBLE
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Datum ellipsoid_in(PG_FUNCTION_ARGS)
#define POLYGONTYPE
Definition: liblwgeom.h:72
double b
Definition: liblwgeom.h:298
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1050
Datum LWGEOM_length2d_ellipsoid(PG_FUNCTION_ARGS)
#define MULTIPOINTTYPE
Definition: liblwgeom.h:73
double radius
Definition: liblwgeom.h:302
void error_if_srid_mismatch(int srid1, int srid2)
Definition: lwutil.c:341
double bigB(double u2)
double lwgeom_length_spheroid(const LWGEOM *geom, const SPHEROID *s)
Calculate the geodetic length of a lwgeom on the unit sphere.
Definition: lwgeodetic.c:2884
Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS)
PG_FUNCTION_INFO_V1(ellipsoid_in)
double mu2(double azimuth, SPHEROID *sphere)
double f
Definition: liblwgeom.h:299
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:61
double distance_ellipse(double lat1, double long1, double lat2, double long2, SPHEROID *sphere)
double distance_sphere_method(double lat1, double long1, double lat2, double long2, SPHEROID *sphere)
double e
Definition: liblwgeom.h:300
char * s
Definition: cu_in_wkt.c:23
double e_sq
Definition: liblwgeom.h:301
Datum distance(PG_FUNCTION_ARGS)
void lwgeom_set_geodetic(LWGEOM *geom, int value)
Set the FLAGS geodetic bit on geometry an all sub-geometries and pointlists.
Definition: lwgeom.c:858
double distance_ellipse_calculation(double lat1, double long1, double lat2, double long2, SPHEROID *sphere)
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:75
double a
Definition: liblwgeom.h:297
#define FALSE
Definition: dbfopen.c:168
char name[20]
Definition: liblwgeom.h:303
double bigA(double u2)
double deltaLongitude(double azimuth, double sigma, double tsm, SPHEROID *sphere)
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:70
Datum geometry_distance_spheroid(PG_FUNCTION_ARGS)
double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance)
Calculate the geodetic distance from lwgeom1 to lwgeom2 on the spheroid.
Definition: lwgeodetic.c:2081
Datum LWGEOM_distance_ellipsoid(PG_FUNCTION_ARGS)
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1297
#define MULTILINETYPE
Definition: liblwgeom.h:74
Datum ellipsoid_out(PG_FUNCTION_ARGS)
#define TRUE
Definition: dbfopen.c:169
int32_t gserialized_get_srid(const GSERIALIZED *s)
Extract the SRID from the serialized form (it is packed into three bytes so this is a handy function)...
Definition: g_serialized.c:69
This library is the generic geometry handling section of PostGIS.