PostGIS  2.1.10dev-r@@SVN_REVISION@@
geography_measurement.c
Go to the documentation of this file.
1 /**********************************************************************
2  * $Id: geography_inout.c 4535 2009-09-28 18:16:21Z colivier $
3  *
4  * PostGIS - Spatial Types for PostgreSQL
5  *
6  * Copyright (C) 2009 Paul Ramsey <pramsey@cleverelephant.ca>
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 #include "../postgis_config.h"
16 
17 #include <math.h>
18 #include <float.h>
19 #include <string.h>
20 #include <stdio.h>
21 #include <errno.h>
22 
23 #include "liblwgeom.h" /* For standard geometry types. */
24 #include "liblwgeom_internal.h" /* For FP comparators. */
25 #include "lwgeom_pg.h" /* For debugging macros. */
26 #include "geography.h" /* For utility functions. */
27 #include "geography_measurement_trees.h" /* For circ_tree caching */
28 #include "lwgeom_transform.h" /* For SRID functions */
29 
30 Datum geography_distance(PG_FUNCTION_ARGS);
31 Datum geography_distance_uncached(PG_FUNCTION_ARGS);
32 Datum geography_distance_tree(PG_FUNCTION_ARGS);
33 Datum geography_dwithin(PG_FUNCTION_ARGS);
34 Datum geography_dwithin_uncached(PG_FUNCTION_ARGS);
35 Datum geography_area(PG_FUNCTION_ARGS);
36 Datum geography_length(PG_FUNCTION_ARGS);
37 Datum geography_expand(PG_FUNCTION_ARGS);
38 Datum geography_point_outside(PG_FUNCTION_ARGS);
39 Datum geography_covers(PG_FUNCTION_ARGS);
40 Datum geography_bestsrid(PG_FUNCTION_ARGS);
41 Datum geography_perimeter(PG_FUNCTION_ARGS);
42 Datum geography_project(PG_FUNCTION_ARGS);
43 Datum geography_azimuth(PG_FUNCTION_ARGS);
44 Datum geography_segmentize(PG_FUNCTION_ARGS);
45 
46 /*
47 ** geography_distance_uncached(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
48 ** returns double distance in meters
49 */
51 Datum geography_distance_uncached(PG_FUNCTION_ARGS)
52 {
53  LWGEOM *lwgeom1 = NULL;
54  LWGEOM *lwgeom2 = NULL;
55  GSERIALIZED *g1 = NULL;
56  GSERIALIZED *g2 = NULL;
57  double distance;
58  /* double tolerance; */
59  bool use_spheroid;
60  SPHEROID s;
61 
62  /* Get our geometry objects loaded into memory. */
63  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
64  g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
65 
66  /* Read our tolerance value. */
67  /* tolerance = PG_GETARG_FLOAT8(2); */
68 
69  /* Read our calculation type. */
70  use_spheroid = PG_GETARG_BOOL(3);
71 
72  /* Initialize spheroid */
73  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g1), &s);
74 
75  /* Set to sphere if requested */
76  if ( ! use_spheroid )
77  s.a = s.b = s.radius;
78 
79  lwgeom1 = lwgeom_from_gserialized(g1);
80  lwgeom2 = lwgeom_from_gserialized(g2);
81 
82  /* Return NULL on empty arguments. */
83  if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
84  {
85  PG_FREE_IF_COPY(g1, 0);
86  PG_FREE_IF_COPY(g2, 1);
87  PG_RETURN_NULL();
88  }
89 
90  /* Make sure we have boxes attached */
91  lwgeom_add_bbox_deep(lwgeom1, NULL);
92  lwgeom_add_bbox_deep(lwgeom2, NULL);
93 
94  distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, FP_TOLERANCE);
95 
96  /* Clean up */
97  lwgeom_free(lwgeom1);
98  lwgeom_free(lwgeom2);
99  PG_FREE_IF_COPY(g1, 0);
100  PG_FREE_IF_COPY(g2, 1);
101 
102  /* Something went wrong, negative return... should already be eloged, return NULL */
103  if ( distance < 0.0 )
104  {
105  PG_RETURN_NULL();
106  }
107 
108  PG_RETURN_FLOAT8(distance);
109 }
110 
111 
112 /*
113 ** geography_distance(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
114 ** returns double distance in meters
115 */
117 Datum geography_distance(PG_FUNCTION_ARGS)
118 {
119  GSERIALIZED* g1 = NULL;
120  GSERIALIZED* g2 = NULL;
121  double distance;
122  double tolerance;
123  bool use_spheroid;
124  SPHEROID s;
125 
126  /* Get our geometry objects loaded into memory. */
127  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
128  g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
129 
130  /* Read our tolerance value. */
131  tolerance = PG_GETARG_FLOAT8(2);
132 
133  /* Read our calculation type. */
134  use_spheroid = PG_GETARG_BOOL(3);
135 
136  /* Initialize spheroid */
137  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g1), &s);
138 
139  /* Set to sphere if requested */
140  if ( ! use_spheroid )
141  s.a = s.b = s.radius;
142 
143  /* Return NULL on empty arguments. */
145  {
146  PG_FREE_IF_COPY(g1, 0);
147  PG_FREE_IF_COPY(g2, 1);
148  PG_RETURN_NULL();
149  }
150 
151  /* Do the brute force calculation if the cached calculation doesn't tick over */
152  if ( LW_FAILURE == geography_distance_cache(fcinfo, g1, g2, &s, &distance) )
153  {
154  LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
155  LWGEOM* lwgeom2 = lwgeom_from_gserialized(g2);
156  distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
157  lwgeom_free(lwgeom1);
158  lwgeom_free(lwgeom2);
159  }
160 
161  /* Clean up */
162  PG_FREE_IF_COPY(g1, 0);
163  PG_FREE_IF_COPY(g2, 1);
164 
165  /* Something went wrong, negative return... should already be eloged, return NULL */
166  if ( distance < 0.0 )
167  {
168  elog(ERROR, "distance returned negative!");
169  PG_RETURN_NULL();
170  }
171 
172  /* Knock off any funny business at the micrometer level, ticket #2168 */
173  distance = round(distance * 10e8) / 10e8;
174 
175  PG_RETURN_FLOAT8(distance);
176 }
177 
178 
179 /*
180 ** geography_dwithin(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
181 ** returns double distance in meters
182 */
184 Datum geography_dwithin(PG_FUNCTION_ARGS)
185 {
186  GSERIALIZED *g1 = NULL;
187  GSERIALIZED *g2 = NULL;
188  double tolerance;
189  double distance;
190  bool use_spheroid;
191  SPHEROID s;
192  int dwithin = LW_FALSE;
193 
194  /* Get our geometry objects loaded into memory. */
195  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
196  g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
197 
198  /* Read our tolerance value. */
199  tolerance = PG_GETARG_FLOAT8(2);
200 
201  /* Read our calculation type. */
202  use_spheroid = PG_GETARG_BOOL(3);
203 
204  /* Initialize spheroid */
205  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g1), &s);
206 
207  /* Set to sphere if requested */
208  if ( ! use_spheroid )
209  s.a = s.b = s.radius;
210 
211  /* Return FALSE on empty arguments. */
213  {
214  PG_FREE_IF_COPY(g1, 0);
215  PG_FREE_IF_COPY(g2, 1);
216  PG_RETURN_BOOL(FALSE);
217  }
218 
219  /* Do the brute force calculation if the cached calculation doesn't tick over */
220  if ( LW_FAILURE == geography_dwithin_cache(fcinfo, g1, g2, &s, tolerance, &dwithin) )
221  {
222  LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
223  LWGEOM* lwgeom2 = lwgeom_from_gserialized(g2);
224  distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
225  /* Something went wrong... */
226  if ( distance < 0.0 )
227  elog(ERROR, "lwgeom_distance_spheroid returned negative!");
228  dwithin = (distance <= tolerance);
229  lwgeom_free(lwgeom1);
230  lwgeom_free(lwgeom2);
231  }
232 
233  /* Clean up */
234  PG_FREE_IF_COPY(g1, 0);
235  PG_FREE_IF_COPY(g2, 1);
236 
237  PG_RETURN_BOOL(dwithin);
238 }
239 
240 
241 /*
242 ** geography_distance_tree(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
243 ** returns double distance in meters
244 */
246 Datum geography_distance_tree(PG_FUNCTION_ARGS)
247 {
248  GSERIALIZED *g1 = NULL;
249  GSERIALIZED *g2 = NULL;
250  double tolerance;
251  double distance;
252  bool use_spheroid;
253  SPHEROID s;
254 
255  /* Get our geometry objects loaded into memory. */
256  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
257  g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
258 
259  /* Return FALSE on empty arguments. */
261  {
262  PG_FREE_IF_COPY(g1, 0);
263  PG_FREE_IF_COPY(g2, 1);
264  PG_RETURN_FLOAT8(0.0);
265  }
266 
267  /* Read our tolerance value. */
268  tolerance = PG_GETARG_FLOAT8(2);
269 
270  /* Read our calculation type. */
271  use_spheroid = PG_GETARG_BOOL(3);
272 
273  /* Initialize spheroid */
274  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g1), &s);
275 
276  /* Set to sphere if requested */
277  if ( ! use_spheroid )
278  s.a = s.b = s.radius;
279 
280  if ( geography_tree_distance(g1, g2, &s, tolerance, &distance) == LW_FAILURE )
281  {
282  elog(ERROR, "geography_distance_tree failed!");
283  PG_RETURN_NULL();
284  }
285 
286  PG_RETURN_FLOAT8(distance);
287 }
288 
289 
290 
291 /*
292 ** geography_dwithin_uncached(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
293 ** returns double distance in meters
294 */
296 Datum geography_dwithin_uncached(PG_FUNCTION_ARGS)
297 {
298  LWGEOM *lwgeom1 = NULL;
299  LWGEOM *lwgeom2 = NULL;
300  GSERIALIZED *g1 = NULL;
301  GSERIALIZED *g2 = NULL;
302  double tolerance;
303  double distance;
304  bool use_spheroid;
305  SPHEROID s;
306 
307  /* Get our geometry objects loaded into memory. */
308  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
309  g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
310 
311  /* Read our tolerance value. */
312  tolerance = PG_GETARG_FLOAT8(2);
313 
314  /* Read our calculation type. */
315  use_spheroid = PG_GETARG_BOOL(3);
316 
317  /* Initialize spheroid */
318  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g1), &s);
319 
320  /* Set to sphere if requested */
321  if ( ! use_spheroid )
322  s.a = s.b = s.radius;
323 
324  lwgeom1 = lwgeom_from_gserialized(g1);
325  lwgeom2 = lwgeom_from_gserialized(g2);
326 
327  /* Return FALSE on empty arguments. */
328  if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
329  {
330  PG_RETURN_BOOL(FALSE);
331  }
332 
333  distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
334 
335  /* Clean up */
336  lwgeom_free(lwgeom1);
337  lwgeom_free(lwgeom2);
338  PG_FREE_IF_COPY(g1, 0);
339  PG_FREE_IF_COPY(g2, 1);
340 
341  /* Something went wrong... should already be eloged, return FALSE */
342  if ( distance < 0.0 )
343  {
344  elog(ERROR, "lwgeom_distance_spheroid returned negative!");
345  PG_RETURN_BOOL(FALSE);
346  }
347 
348  PG_RETURN_BOOL(distance <= tolerance);
349 }
350 
351 
352 /*
353 ** geography_expand(GSERIALIZED *g) returns *GSERIALIZED
354 **
355 ** warning, this tricky little function does not expand the
356 ** geometry at all, just re-writes bounding box value to be
357 ** a bit bigger. only useful when passing the result along to
358 ** an index operator (&&)
359 */
361 Datum geography_expand(PG_FUNCTION_ARGS)
362 {
363  GSERIALIZED *g = NULL;
364  GSERIALIZED *g_out = NULL;
365  double distance;
366 
367  /* Get a wholly-owned pointer to the geography */
368  g = (GSERIALIZED*)PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
369 
370  /* Read our distance value and normalize to unit-sphere. */
371  distance = PG_GETARG_FLOAT8(1) / WGS84_RADIUS;
372 
373  /* Try the expansion */
374  g_out = gserialized_expand(g, distance);
375 
376  /* If the expansion fails, the return our input */
377  if ( g_out == NULL )
378  {
379  PG_RETURN_POINTER(g);
380  }
381 
382  if ( g_out != g )
383  {
384  pfree(g);
385  }
386 
387  PG_RETURN_POINTER(g_out);
388 }
389 
390 /*
391 ** geography_area(GSERIALIZED *g)
392 ** returns double area in meters square
393 */
395 Datum geography_area(PG_FUNCTION_ARGS)
396 {
397  LWGEOM *lwgeom = NULL;
398  GSERIALIZED *g = NULL;
399  GBOX gbox;
400  double area;
401  bool use_spheroid = LW_TRUE;
402  SPHEROID s;
403 
404  /* Get our geometry object loaded into memory. */
405  g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
406 
407  /* Read our calculation type */
408  use_spheroid = PG_GETARG_BOOL(1);
409 
410  /* Initialize spheroid */
411  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g), &s);
412 
413  lwgeom = lwgeom_from_gserialized(g);
414 
415  /* EMPTY things have no area */
416  if ( lwgeom_is_empty(lwgeom) )
417  {
418  lwgeom_free(lwgeom);
419  PG_RETURN_FLOAT8(0.0);
420  }
421 
422  if ( lwgeom->bbox )
423  gbox = *(lwgeom->bbox);
424  else
425  lwgeom_calculate_gbox_geodetic(lwgeom, &gbox);
426 
427  /* Test for cases that are currently not handled by spheroid code */
428  if ( use_spheroid )
429  {
430  /* We can't circle the poles right now */
431  if ( FP_GTEQ(gbox.zmax,1.0) || FP_LTEQ(gbox.zmin,-1.0) )
432  use_spheroid = LW_FALSE;
433  /* We can't cross the equator right now */
434  if ( gbox.zmax > 0.0 && gbox.zmin < 0.0 )
435  use_spheroid = LW_FALSE;
436  }
437 
438  /* User requests spherical calculation, turn our spheroid into a sphere */
439  if ( ! use_spheroid )
440  s.a = s.b = s.radius;
441 
442  /* Calculate the area */
443  if ( use_spheroid )
444  area = lwgeom_area_spheroid(lwgeom, &s);
445  else
446  area = lwgeom_area_sphere(lwgeom, &s);
447 
448  /* Clean up */
449  lwgeom_free(lwgeom);
450  PG_FREE_IF_COPY(g, 0);
451 
452  /* Something went wrong... */
453  if ( area < 0.0 )
454  {
455  elog(ERROR, "lwgeom_area_spher(oid) returned area < 0.0");
456  PG_RETURN_NULL();
457  }
458 
459  PG_RETURN_FLOAT8(area);
460 }
461 
462 /*
463 ** geography_perimeter(GSERIALIZED *g)
464 ** returns double perimeter in meters for area features
465 */
467 Datum geography_perimeter(PG_FUNCTION_ARGS)
468 {
469  LWGEOM *lwgeom = NULL;
470  GSERIALIZED *g = NULL;
471  double length;
472  bool use_spheroid = LW_TRUE;
473  SPHEROID s;
474  int type;
475 
476  /* Get our geometry object loaded into memory. */
477  g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
478 
479  /* Only return for area features. */
480  type = gserialized_get_type(g);
481  if ( ! (type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE) )
482  {
483  PG_RETURN_FLOAT8(0.0);
484  }
485 
486  lwgeom = lwgeom_from_gserialized(g);
487 
488  /* EMPTY things have no perimeter */
489  if ( lwgeom_is_empty(lwgeom) )
490  {
491  lwgeom_free(lwgeom);
492  PG_RETURN_FLOAT8(0.0);
493  }
494 
495  /* Read our calculation type */
496  use_spheroid = PG_GETARG_BOOL(1);
497 
498  /* Initialize spheroid */
499  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g), &s);
500 
501  /* User requests spherical calculation, turn our spheroid into a sphere */
502  if ( ! use_spheroid )
503  s.a = s.b = s.radius;
504 
505  /* Calculate the length */
506  length = lwgeom_length_spheroid(lwgeom, &s);
507 
508  /* Something went wrong... */
509  if ( length < 0.0 )
510  {
511  elog(ERROR, "lwgeom_length_spheroid returned length < 0.0");
512  PG_RETURN_NULL();
513  }
514 
515  /* Clean up, but not all the way to the point arrays */
516  lwgeom_free(lwgeom);
517 
518  PG_FREE_IF_COPY(g, 0);
519  PG_RETURN_FLOAT8(length);
520 }
521 
522 /*
523 ** geography_length(GSERIALIZED *g)
524 ** returns double length in meters
525 */
527 Datum geography_length(PG_FUNCTION_ARGS)
528 {
529  LWGEOM *lwgeom = NULL;
530  GSERIALIZED *g = NULL;
531  double length;
532  bool use_spheroid = LW_TRUE;
533  SPHEROID s;
534 
535  /* Get our geometry object loaded into memory. */
536  g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
537  lwgeom = lwgeom_from_gserialized(g);
538 
539  /* EMPTY things have no length */
540  if ( lwgeom_is_empty(lwgeom) || lwgeom->type == POLYGONTYPE || lwgeom->type == MULTIPOLYGONTYPE )
541  {
542  lwgeom_free(lwgeom);
543  PG_RETURN_FLOAT8(0.0);
544  }
545 
546  /* Read our calculation type */
547  use_spheroid = PG_GETARG_BOOL(1);
548 
549  /* Initialize spheroid */
550  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g), &s);
551 
552  /* User requests spherical calculation, turn our spheroid into a sphere */
553  if ( ! use_spheroid )
554  s.a = s.b = s.radius;
555 
556  /* Calculate the length */
557  length = lwgeom_length_spheroid(lwgeom, &s);
558 
559  /* Something went wrong... */
560  if ( length < 0.0 )
561  {
562  elog(ERROR, "lwgeom_length_spheroid returned length < 0.0");
563  PG_RETURN_NULL();
564  }
565 
566  /* Clean up */
567  lwgeom_free(lwgeom);
568 
569  PG_FREE_IF_COPY(g, 0);
570  PG_RETURN_FLOAT8(length);
571 }
572 
573 /*
574 ** geography_point_outside(GSERIALIZED *g)
575 ** returns point outside the object
576 */
578 Datum geography_point_outside(PG_FUNCTION_ARGS)
579 {
580  GBOX gbox;
581  GSERIALIZED *g = NULL;
582  GSERIALIZED *g_out = NULL;
583  size_t g_out_size;
584  LWPOINT *lwpoint = NULL;
585  POINT2D pt;
586 
587  /* Get our geometry object loaded into memory. */
588  g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
589 
590  /* We need the bounding box to get an outside point for area algorithm */
591  if ( gserialized_get_gbox_p(g, &gbox) == LW_FAILURE )
592  {
593  POSTGIS_DEBUG(4,"gserialized_get_gbox_p returned LW_FAILURE");
594  elog(ERROR, "Error in gserialized_get_gbox_p calculation.");
595  PG_RETURN_NULL();
596  }
597  POSTGIS_DEBUGF(4, "got gbox %s", gbox_to_string(&gbox));
598 
599  /* Get an exterior point, based on this gbox */
600  gbox_pt_outside(&gbox, &pt);
601 
602  lwpoint = lwpoint_make2d(4326, pt.x, pt.y);
603 
604  g_out = gserialized_from_lwgeom((LWGEOM*)lwpoint, 1, &g_out_size);
605  SET_VARSIZE(g_out, g_out_size);
606 
607  PG_FREE_IF_COPY(g, 0);
608  PG_RETURN_POINTER(g_out);
609 
610 }
611 
612 /*
613 ** geography_covers(GSERIALIZED *g, GSERIALIZED *g) returns boolean
614 ** Only works for (multi)points and (multi)polygons currently.
615 ** Attempts a simple point-in-polygon test on the polygon and point.
616 ** Current algorithm does not distinguish between points on edge
617 ** and points within.
618 */
620 Datum geography_covers(PG_FUNCTION_ARGS)
621 {
622  LWGEOM *lwgeom1 = NULL;
623  LWGEOM *lwgeom2 = NULL;
624  GSERIALIZED *g1 = NULL;
625  GSERIALIZED *g2 = NULL;
626  int type1, type2;
627  int result = LW_FALSE;
628 
629  /* Get our geometry objects loaded into memory. */
630  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
631  g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
632 
633  type1 = gserialized_get_type(g1);
634  type2 = gserialized_get_type(g2);
635 
636  /* Right now we only handle points and polygons */
637  if ( ! ( (type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE || type1 == COLLECTIONTYPE) &&
638  (type2 == POINTTYPE || type2 == MULTIPOINTTYPE || type2 == COLLECTIONTYPE) ) )
639  {
640  elog(ERROR, "geography_covers: only POLYGON and POINT types are currently supported");
641  PG_RETURN_NULL();
642  }
643 
644  /* Construct our working geometries */
645  lwgeom1 = lwgeom_from_gserialized(g1);
646  lwgeom2 = lwgeom_from_gserialized(g2);
647 
648  /* EMPTY never intersects with another geometry */
649  if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
650  {
651  lwgeom_free(lwgeom1);
652  lwgeom_free(lwgeom2);
653  PG_FREE_IF_COPY(g1, 0);
654  PG_FREE_IF_COPY(g2, 1);
655  PG_RETURN_BOOL(false);
656  }
657 
658  /* Calculate answer */
659  result = lwgeom_covers_lwgeom_sphere(lwgeom1, lwgeom2);
660 
661  /* Clean up */
662  lwgeom_free(lwgeom1);
663  lwgeom_free(lwgeom2);
664  PG_FREE_IF_COPY(g1, 0);
665  PG_FREE_IF_COPY(g2, 1);
666 
667  PG_RETURN_BOOL(result);
668 }
669 
670 /*
671 ** geography_bestsrid(GSERIALIZED *g, GSERIALIZED *g) returns int
672 ** Utility function. Returns negative SRID numbers that match to the
673 ** numbers handled in code by the transform(lwgeom, srid) function.
674 ** UTM, polar stereographic and mercator as fallback. To be used
675 ** in wrapping existing geometry functions in SQL to provide access
676 ** to them in the geography module.
677 */
679 Datum geography_bestsrid(PG_FUNCTION_ARGS)
680 {
681  GBOX gbox, gbox1, gbox2;
682  GSERIALIZED *g1 = NULL;
683  GSERIALIZED *g2 = NULL;
684  int empty1 = LW_FALSE;
685  int empty2 = LW_FALSE;
686  double xwidth, ywidth;
687  POINT2D center;
688 
689  Datum d1 = PG_GETARG_DATUM(0);
690  Datum d2 = PG_GETARG_DATUM(1);
691 
692  /* Get our geometry objects loaded into memory. */
693  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(d1);
694  /* Synchronize our box types */
695  gbox1.flags = g1->flags;
696  /* Calculate if the geometry is empty. */
697  empty1 = gserialized_is_empty(g1);
698  /* Calculate a geocentric bounds for the objects */
699  if ( ! empty1 && gserialized_get_gbox_p(g1, &gbox1) == LW_FAILURE )
700  elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g1, &gbox1)");
701 
702  POSTGIS_DEBUGF(4, "calculated gbox = %s", gbox_to_string(&gbox1));
703 
704  /* If we have a unique second argument, fill in all the necessary variables. */
705  if ( d1 != d2 )
706  {
707  g2 = (GSERIALIZED*)PG_DETOAST_DATUM(d2);
708  gbox2.flags = g2->flags;
709  empty2 = gserialized_is_empty(g2);
710  if ( ! empty2 && gserialized_get_gbox_p(g2, &gbox2) == LW_FAILURE )
711  elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g2, &gbox2)");
712  }
713  /*
714  ** If no unique second argument, copying the box for the first
715  ** argument will give us the right answer for all subsequent tests.
716  */
717  else
718  {
719  gbox = gbox2 = gbox1;
720  }
721 
722  /* Both empty? We don't have an answer. */
723  if ( empty1 && empty2 )
724  PG_RETURN_NULL();
725 
726  /* One empty? We can use the other argument values as infill. Otherwise merge the boxen */
727  if ( empty1 )
728  gbox = gbox2;
729  else if ( empty2 )
730  gbox = gbox1;
731  else
732  gbox_union(&gbox1, &gbox2, &gbox);
733 
734  gbox_centroid(&gbox, &center);
735 
736  /* Width and height in degrees */
737  xwidth = 180.0 * gbox_angular_width(&gbox) / M_PI;
738  ywidth = 180.0 * gbox_angular_height(&gbox) / M_PI;
739 
740  POSTGIS_DEBUGF(2, "xwidth %g", xwidth);
741  POSTGIS_DEBUGF(2, "ywidth %g", ywidth);
742  POSTGIS_DEBUGF(2, "center POINT(%g %g)", center.x, center.y);
743 
744  /* Are these data arctic? Lambert Azimuthal Equal Area North. */
745  if ( center.y > 70.0 && ywidth < 45.0 )
746  {
747  PG_RETURN_INT32(SRID_NORTH_LAMBERT);
748  }
749 
750  /* Are these data antarctic? Lambert Azimuthal Equal Area South. */
751  if ( center.y < -70.0 && ywidth < 45.0 )
752  {
753  PG_RETURN_INT32(SRID_SOUTH_LAMBERT);
754  }
755 
756  /*
757  ** Can we fit these data into one UTM zone?
758  ** We will assume we can push things as
759  ** far as a half zone past a zone boundary.
760  ** Note we have no handling for the date line in here.
761  */
762  if ( xwidth < 6.0 )
763  {
764  int zone = floor((center.x + 180.0) / 6.0);
765 
766  if ( zone > 59 ) zone = 59;
767 
768  /* Are these data below the equator? UTM South. */
769  if ( center.y < 0.0 )
770  {
771  PG_RETURN_INT32( SRID_SOUTH_UTM_START + zone );
772  }
773  /* Are these data above the equator? UTM North. */
774  else
775  {
776  PG_RETURN_INT32( SRID_NORTH_UTM_START + zone );
777  }
778  }
779 
780  /*
781  ** Can we fit into a custom LAEA area? (30 degrees high, variable width)
782  ** We will allow overlap into adjoining areas, but use a slightly narrower test (25) to try
783  ** and minimize the worst case.
784  ** Again, we are hoping the dateline doesn't trip us up much
785  */
786  if ( ywidth < 25.0 )
787  {
788  int xzone = -1;
789  int yzone = 3 + floor(center.y / 30.0); /* (range of 0-5) */
790 
791  /* Equatorial band, 12 zones, 30 degrees wide */
792  if ( (yzone == 2 || yzone == 3) && xwidth < 30.0 )
793  {
794  xzone = 6 + floor(center.x / 30.0);
795  }
796  /* Temperate band, 8 zones, 45 degrees wide */
797  else if ( (yzone == 1 || yzone == 4) && xwidth < 45.0 )
798  {
799  xzone = 4 + floor(center.x / 45.0);
800  }
801  /* Arctic band, 4 zones, 90 degrees wide */
802  else if ( (yzone == 0 || yzone == 5) && xwidth < 90.0 )
803  {
804  xzone = 2 + floor(center.x / 90.0);
805  }
806 
807  /* Did we fit into an appropriate xzone? */
808  if ( xzone != -1 )
809  {
810  PG_RETURN_INT32(SRID_LAEA_START + 20 * yzone + xzone);
811  }
812  }
813 
814  /*
815  ** Running out of options... fall-back to Mercator
816  ** and hope for the best.
817  */
818  PG_RETURN_INT32(SRID_WORLD_MERCATOR);
819 
820 }
821 
822 /*
823 ** geography_project(GSERIALIZED *g, distance, azimuth)
824 ** returns point of projection given start point,
825 ** azimuth in radians (bearing) and distance in meters
826 */
828 Datum geography_project(PG_FUNCTION_ARGS)
829 {
830  LWGEOM *lwgeom = NULL;
831  LWPOINT *lwp_projected;
832  GSERIALIZED *g = NULL;
833  GSERIALIZED *g_out = NULL;
834  double azimuth, distance;
835  SPHEROID s;
836  uint32_t type;
837 
838  /* Return NULL on NULL distance or geography */
839  if ( PG_ARGISNULL(0) || PG_ARGISNULL(1) )
840  PG_RETURN_NULL();
841 
842  /* Get our geometry object loaded into memory. */
843  g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
844 
845  /* Only return for points. */
846  type = gserialized_get_type(g);
847  if ( type != POINTTYPE )
848  {
849  elog(ERROR, "ST_Project(geography) is only valid for point inputs");
850  PG_RETURN_NULL();
851  }
852 
853  distance = PG_GETARG_FLOAT8(1); /* Distance in Meters */
854  lwgeom = lwgeom_from_gserialized(g);
855 
856  /* EMPTY things cannot be projected from */
857  if ( lwgeom_is_empty(lwgeom) )
858  {
859  lwgeom_free(lwgeom);
860  elog(ERROR, "ST_Project(geography) cannot project from an empty start point");
861  PG_RETURN_NULL();
862  }
863 
864  if ( PG_ARGISNULL(2) )
865  azimuth = 0.0;
866  else
867  azimuth = PG_GETARG_FLOAT8(2); /* Azimuth in Radians */
868 
869  /* Initialize spheroid */
870  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g), &s);
871 
872  /* Handle the zero distance case */
873  if( FP_EQUALS(distance, 0.0) )
874  {
875  PG_RETURN_POINTER(g);
876  }
877 
878  /* Calculate the length */
879  lwp_projected = lwgeom_project_spheroid(lwgeom_as_lwpoint(lwgeom), &s, distance, azimuth);
880 
881  /* Something went wrong... */
882  if ( lwp_projected == NULL )
883  {
884  elog(ERROR, "lwgeom_project_spheroid returned null");
885  PG_RETURN_NULL();
886  }
887 
888  /* Clean up, but not all the way to the point arrays */
889  lwgeom_free(lwgeom);
890  g_out = geography_serialize(lwpoint_as_lwgeom(lwp_projected));
891  lwpoint_free(lwp_projected);
892 
893  PG_FREE_IF_COPY(g, 0);
894  PG_RETURN_POINTER(g_out);
895 }
896 
897 
898 /*
899 ** geography_azimuth(GSERIALIZED *g1, GSERIALIZED *g2)
900 ** returns direction between points (north = 0)
901 ** azimuth (bearing) and distance
902 */
904 Datum geography_azimuth(PG_FUNCTION_ARGS)
905 {
906  LWGEOM *lwgeom1 = NULL;
907  LWGEOM *lwgeom2 = NULL;
908  GSERIALIZED *g1 = NULL;
909  GSERIALIZED *g2 = NULL;
910  double azimuth;
911  SPHEROID s;
912  uint32_t type1, type2;
913 
914  /* Get our geometry object loaded into memory. */
915  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
916  g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
917 
918  /* Only return for points. */
919  type1 = gserialized_get_type(g1);
920  type2 = gserialized_get_type(g2);
921  if ( type1 != POINTTYPE || type2 != POINTTYPE )
922  {
923  elog(ERROR, "ST_Azimuth(geography, geography) is only valid for point inputs");
924  PG_RETURN_NULL();
925  }
926 
927  lwgeom1 = lwgeom_from_gserialized(g1);
928  lwgeom2 = lwgeom_from_gserialized(g2);
929 
930  /* EMPTY things cannot be used */
931  if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
932  {
933  lwgeom_free(lwgeom1);
934  lwgeom_free(lwgeom2);
935  elog(ERROR, "ST_Azimuth(geography, geography) cannot work with empty points");
936  PG_RETURN_NULL();
937  }
938 
939  /* Initialize spheroid */
940  spheroid_init_from_srid(fcinfo, gserialized_get_srid(g1), &s);
941 
942  /* Calculate the direction */
943  azimuth = lwgeom_azumith_spheroid(lwgeom_as_lwpoint(lwgeom1), lwgeom_as_lwpoint(lwgeom2), &s);
944 
945  /* Clean up */
946  lwgeom_free(lwgeom1);
947  lwgeom_free(lwgeom2);
948 
949  PG_FREE_IF_COPY(g1, 0);
950  PG_FREE_IF_COPY(g2, 1);
951 
952  /* Return NULL for unknown (same point) azimuth */
953  if( isnan(azimuth) )
954  {
955  PG_RETURN_NULL();
956  }
957 
958  PG_RETURN_FLOAT8(azimuth);
959 }
960 
961 
962 
963 /*
964 ** geography_segmentize(GSERIALIZED *g1, double max_seg_length)
965 ** returns densified geometry with no segment longer than max
966 */
968 Datum geography_segmentize(PG_FUNCTION_ARGS)
969 {
970  LWGEOM *lwgeom1 = NULL;
971  LWGEOM *lwgeom2 = NULL;
972  GSERIALIZED *g1 = NULL;
973  GSERIALIZED *g2 = NULL;
974  double max_seg_length;
975  uint32_t type1;
976 
977  /* Get our geometry object loaded into memory. */
978  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
979  type1 = gserialized_get_type(g1);
980 
981  /* Convert max_seg_length from metric units to radians */
982  max_seg_length = PG_GETARG_FLOAT8(1) / WGS84_RADIUS;
983 
984  /* We can't densify points or points, reflect them back */
985  if ( type1 == POINTTYPE || type1 == MULTIPOINTTYPE || gserialized_is_empty(g1) )
986  PG_RETURN_POINTER(g1);
987 
988  /* Deserialize */
989  lwgeom1 = lwgeom_from_gserialized(g1);
990 
991  /* Calculate the densified geometry */
992  lwgeom2 = lwgeom_segmentize_sphere(lwgeom1, max_seg_length);
993 
994  /*
995  ** Set the geodetic flag so subsequent
996  ** functions do the right thing.
997  */
998  lwgeom_set_geodetic(lwgeom2, true);
999 
1000  /* Recalculate the boxes after re-setting the geodetic bit */
1001  lwgeom_drop_bbox(lwgeom2);
1002  lwgeom_add_bbox(lwgeom2);
1003 
1004  g2 = geography_serialize(lwgeom2);
1005 
1006  /* Clean up */
1007  lwgeom_free(lwgeom1);
1008  lwgeom_free(lwgeom2);
1009  PG_FREE_IF_COPY(g1, 0);
1010 
1011  PG_RETURN_POINTER(g2);
1012 }
1013 
1014 
int gserialized_get_gbox_p(const GSERIALIZED *g, GBOX *box)
Read the bounding box off a serialization and calculate one if it is not already there.
Definition: g_serialized.c:373
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:56
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:2612
void lwgeom_add_bbox_deep(LWGEOM *lwgeom, GBOX *gbox)
Compute a box for geom and all sub-geometries, if not already computed.
Definition: lwgeom.c:567
GBOX * bbox
Definition: liblwgeom.h:354
Datum geography_segmentize(PG_FUNCTION_ARGS)
double lwgeom_azumith_spheroid(const LWPOINT *r, const LWPOINT *s, const SPHEROID *spheroid)
Calculate the bearing between two points on a spheroid.
Definition: lwgeodetic.c:2047
Datum geography_distance(PG_FUNCTION_ARGS)
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
char * gbox_to_string(const GBOX *gbox)
Allocate a string representation of the GBOX, based on dimensionality of flags.
Definition: g_box.c:328
#define WGS84_RADIUS
Definition: liblwgeom.h:98
Datum geography_length(PG_FUNCTION_ARGS)
int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
Calculate covers predicate for two lwgeoms on the sphere.
Definition: lwgeodetic.c:2300
#define POLYGONTYPE
Definition: liblwgeom.h:62
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:130
Datum area(PG_FUNCTION_ARGS)
double b
Definition: liblwgeom.h:270
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:180
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1006
Datum geography_dwithin_uncached(PG_FUNCTION_ARGS)
#define MULTIPOINTTYPE
Definition: liblwgeom.h:63
double radius
Definition: liblwgeom.h:274
int gbox_centroid(const GBOX *gbox, POINT2D *out)
Computes the average(ish) center of the box and returns success.
Definition: lwgeodetic.c:244
#define FP_LTEQ(A, B)
LWPOINT * lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, double distance, double azimuth)
Calculate the location of a point on a spheroid, give a start point, bearing and distance.
Definition: lwgeodetic.c:1992
int geography_tree_distance(const GSERIALIZED *g1, const GSERIALIZED *g2, const SPHEROID *s, double tolerance, double *distance)
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:80
double lwgeom_length_spheroid(const LWGEOM *geom, const SPHEROID *s)
Calculate the geodetic length of a lwgeom on the unit sphere.
Definition: lwgeodetic.c:2881
Datum geography_area(PG_FUNCTION_ARGS)
char ** result
Definition: liblwgeom.h:218
void lwgeom_drop_bbox(LWGEOM *lwgeom)
Call this function to drop BBOX and SRID from LWGEOM.
Definition: lwgeom.c:542
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: g_serialized.c:140
#define LW_FAILURE
Definition: liblwgeom.h:54
GSERIALIZED * gserialized_expand(GSERIALIZED *g, double distance)
Return a GSERIALIZED with an expanded bounding box.
double x
Definition: liblwgeom.h:284
Datum geography_bestsrid(PG_FUNCTION_ARGS)
double zmax
Definition: liblwgeom.h:253
#define LW_FALSE
Definition: liblwgeom.h:52
PG_FUNCTION_INFO_V1(geography_distance_uncached)
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:51
Datum geography_point_outside(PG_FUNCTION_ARGS)
GSERIALIZED * gserialized_from_lwgeom(LWGEOM *geom, int is_geodetic, size_t *size)
Allocate a new GSERIALIZED from an LWGEOM.
Definition: g_serialized.c:908
#define FP_TOLERANCE
Floating point comparitors.
double y
Definition: liblwgeom.h:284
char * s
Definition: cu_in_wkt.c:24
Datum distance(PG_FUNCTION_ARGS)
int gbox_union(const GBOX *g1, const GBOX *g2, GBOX *gout)
Update the output GBOX to be large enough to include both inputs.
Definition: g_box.c:109
double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the geodetic area of a lwgeom on the sphere.
Definition: lwgeodetic.c:1924
void lwgeom_set_geodetic(LWGEOM *geom, int value)
Set the FLAGS geodetic bit on geometry an all sub-geometries and pointlists.
Definition: lwgeom.c:814
uint8_t flags
Definition: liblwgeom.h:247
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:65
void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
Calculate a spherical point that falls outside the geocentric gbox.
Definition: lwgeodetic.c:1443
double a
Definition: liblwgeom.h:269
int geography_dwithin_cache(FunctionCallInfoData *fcinfo, const GSERIALIZED *g1, const GSERIALIZED *g2, const SPHEROID *s, double tolerance, int *dwithin)
#define FP_GTEQ(A, B)
LWGEOM * lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length)
Derive a new geometry with vertices added to ensure no vertex is more than max_seg_length (in radians...
Definition: lwgeodetic.c:1638
#define FALSE
Definition: dbfopen.c:169
Datum geography_expand(PG_FUNCTION_ARGS)
Datum geography_distance_tree(PG_FUNCTION_ARGS)
double zmin
Definition: liblwgeom.h:252
double gbox_angular_height(const GBOX *gbox)
GBOX utility functions to figure out coverage/location on the globe.
Definition: lwgeodetic.c:165
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:555
int geography_distance_cache(FunctionCallInfoData *fcinfo, const GSERIALIZED *g1, const GSERIALIZED *g2, const SPHEROID *s, double *distance)
Datum geography_covers(PG_FUNCTION_ARGS)
Datum geography_azimuth(PG_FUNCTION_ARGS)
uint8_t type
Definition: liblwgeom.h:352
#define FP_EQUALS(A, B)
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:254
Datum geography_distance_uncached(PG_FUNCTION_ARGS)
double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the geodetic area of a lwgeom on the spheroid.
Definition: lwspheroid.c:513
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:2078
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:1229
Datum geography_dwithin(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:70
#define COLLECTIONTYPE
Definition: liblwgeom.h:66
This library is the generic geometry handling section of PostGIS.
uint8_t flags
Definition: liblwgeom.h:339
Datum geography_perimeter(PG_FUNCTION_ARGS)
Datum geography_project(PG_FUNCTION_ARGS)
double gbox_angular_width(const GBOX *gbox)
Returns the angular width (longitudinal span) of the box in radians.
Definition: lwgeodetic.c:192