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