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