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