PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches
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
51Datum geography_distance(PG_FUNCTION_ARGS);
52Datum geography_distance_uncached(PG_FUNCTION_ARGS);
53Datum geography_distance_knn(PG_FUNCTION_ARGS);
54Datum geography_distance_tree(PG_FUNCTION_ARGS);
55Datum geography_dwithin(PG_FUNCTION_ARGS);
56Datum geography_dwithin_uncached(PG_FUNCTION_ARGS);
57Datum geography_area(PG_FUNCTION_ARGS);
58Datum geography_length(PG_FUNCTION_ARGS);
59Datum geography_expand(PG_FUNCTION_ARGS);
60Datum geography_point_outside(PG_FUNCTION_ARGS);
61Datum geography_covers(PG_FUNCTION_ARGS);
62Datum geography_coveredby(PG_FUNCTION_ARGS);
63Datum geography_bestsrid(PG_FUNCTION_ARGS);
64Datum geography_perimeter(PG_FUNCTION_ARGS);
65Datum geography_project(PG_FUNCTION_ARGS);
66Datum geography_azimuth(PG_FUNCTION_ARGS);
67Datum geography_segmentize(PG_FUNCTION_ARGS);
68
69Datum geography_line_locate_point(PG_FUNCTION_ARGS);
70Datum geography_line_interpolate_point(PG_FUNCTION_ARGS);
71Datum geography_line_substring(PG_FUNCTION_ARGS);
72Datum geography_closestpoint(PG_FUNCTION_ARGS);
73Datum geography_shortestline(PG_FUNCTION_ARGS);
74
75
76PG_FUNCTION_INFO_V1(geography_distance_knn);
77Datum geography_distance_knn(PG_FUNCTION_ARGS)
78{
79 LWGEOM *lwgeom1 = NULL;
80 LWGEOM *lwgeom2 = NULL;
81 GSERIALIZED *g1 = NULL;
82 GSERIALIZED *g2 = NULL;
83 double distance;
84 double tolerance = FP_TOLERANCE;
85 bool use_spheroid = false; /* must use sphere, can't get index to harmonize with spheroid */
86 SPHEROID s;
87
88 /* Get our geometry objects loaded into memory. */
89 g1 = PG_GETARG_GSERIALIZED_P(0);
90 g2 = PG_GETARG_GSERIALIZED_P(1);
91
92 gserialized_error_if_srid_mismatch(g1, g2, __func__);
93
94 /* Initialize spheroid */
95 spheroid_init_from_srid(gserialized_get_srid(g1), &s);
96
97 /* Set to sphere if requested */
98 if ( ! use_spheroid )
99 s.a = s.b = s.radius;
100
101 lwgeom1 = lwgeom_from_gserialized(g1);
102 lwgeom2 = lwgeom_from_gserialized(g2);
103
104 /* Return NULL on empty arguments. */
105 if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
106 {
107 PG_FREE_IF_COPY(g1, 0);
108 PG_FREE_IF_COPY(g2, 1);
109 PG_RETURN_NULL();
110 }
111
112 /* Make sure we have boxes attached */
113 lwgeom_add_bbox_deep(lwgeom1, NULL);
114 lwgeom_add_bbox_deep(lwgeom2, NULL);
115
116 distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
117
118 POSTGIS_DEBUGF(2, "[GIST] '%s' got distance %g", __func__, distance);
119
120 /* Clean up */
121 lwgeom_free(lwgeom1);
122 lwgeom_free(lwgeom2);
123 PG_FREE_IF_COPY(g1, 0);
124 PG_FREE_IF_COPY(g2, 1);
125
126 /* Something went wrong, negative return... should already be eloged, return NULL */
127 if ( distance < 0.0 )
128 {
129 PG_RETURN_NULL();
130 }
131
132 PG_RETURN_FLOAT8(distance);
133}
134
135/*
136** geography_distance_uncached(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
137** returns double distance in meters
138*/
139PG_FUNCTION_INFO_V1(geography_distance_uncached);
140Datum geography_distance_uncached(PG_FUNCTION_ARGS)
141{
142 LWGEOM *lwgeom1 = NULL;
143 LWGEOM *lwgeom2 = NULL;
144 GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
145 GSERIALIZED *g2 = PG_GETARG_GSERIALIZED_P(1);
146 double distance;
147 double tolerance = FP_TOLERANCE;
148 bool use_spheroid = true;
149 SPHEROID s;
150
151 /* Read our tolerance value. */
152 if ( PG_NARGS() > 2 && ! PG_ARGISNULL(2) )
153 tolerance = PG_GETARG_FLOAT8(2);
154
155 /* Read our calculation type. */
156 if ( PG_NARGS() > 3 && ! PG_ARGISNULL(3) )
157 use_spheroid = PG_GETARG_BOOL(3);
158
159 gserialized_error_if_srid_mismatch(g1, g2, __func__);
160
161 /* Initialize spheroid */
162 spheroid_init_from_srid(gserialized_get_srid(g1), &s);
163
164 /* Set to sphere if requested */
165 if ( ! use_spheroid )
166 s.a = s.b = s.radius;
167
168 lwgeom1 = lwgeom_from_gserialized(g1);
169 lwgeom2 = lwgeom_from_gserialized(g2);
170
171 /* Return NULL on empty arguments. */
172 if ( !lwgeom1 || !lwgeom2 || lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
173 {
174 PG_FREE_IF_COPY(g1, 0);
175 PG_FREE_IF_COPY(g2, 1);
176 PG_RETURN_NULL();
177 }
178
179 /* Make sure we have boxes attached */
180 lwgeom_add_bbox_deep(lwgeom1, NULL);
181 lwgeom_add_bbox_deep(lwgeom2, NULL);
182
183 distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
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*/
205PG_FUNCTION_INFO_V1(geography_distance);
206Datum 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 if (PG_NARGS() > 2)
217 use_spheroid = PG_GETARG_BOOL(2);
218
219 gserialized_error_if_srid_mismatch(g1, g2, __func__);
220
221 /* Initialize spheroid */
222 spheroid_init_from_srid(gserialized_get_srid(g1), &s);
223
224 /* Set to sphere if requested */
225 if ( ! use_spheroid )
226 s.a = s.b = s.radius;
227
228 /* Return NULL on empty arguments. */
229 if ( gserialized_is_empty(g1) || gserialized_is_empty(g2) )
230 {
231 PG_RETURN_NULL();
232 }
233
234 /* Do the brute force calculation if the cached calculation doesn't tick over */
235 if (LW_FAILURE == geography_distance_cache(fcinfo, shared_geom1, shared_geom2, &s, &distance))
236 {
237 /* default to using tree-based distance calculation at all times */
238 /* in standard distance call. */
239 geography_tree_distance(g1, g2, &s, FP_TOLERANCE, &distance);
240 /*
241 LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
242 LWGEOM* lwgeom2 = lwgeom_from_gserialized(g2);
243 distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
244 lwgeom_free(lwgeom1);
245 lwgeom_free(lwgeom2);
246 */
247 }
248
249 /* Knock off any funny business at the nanometer level, ticket #2168 */
250 distance = round(distance * INVMINDIST) / INVMINDIST;
251
252 /* Something went wrong, negative return... should already be eloged, return NULL */
253 if ( distance < 0.0 )
254 {
255 elog(ERROR, "distance returned negative!");
256 PG_RETURN_NULL();
257 }
258
259 PG_RETURN_FLOAT8(distance);
260}
261
262/*
263** geography_dwithin(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
264** returns double distance in meters
265*/
266PG_FUNCTION_INFO_V1(geography_dwithin);
267Datum geography_dwithin(PG_FUNCTION_ARGS)
268{
269 SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
270 SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
271 const GSERIALIZED *g1 = shared_gserialized_get(shared_geom1);
272 const GSERIALIZED *g2 = shared_gserialized_get(shared_geom2);
273 SPHEROID s;
274 double tolerance = FP_TOLERANCE;
275 bool use_spheroid = true;
276 double distance;
277 int dwithin = LW_FALSE;
278
279 gserialized_error_if_srid_mismatch(g1, g2, __func__);
280
281 /* Read our tolerance value. */
282 if ( PG_NARGS() > 2 && ! PG_ARGISNULL(2) )
283 tolerance = PG_GETARG_FLOAT8(2);
284
285 /* Read our calculation type. */
286 if ( PG_NARGS() > 3 && ! PG_ARGISNULL(3) )
287 use_spheroid = PG_GETARG_BOOL(3);
288
289 /* Initialize spheroid */
290 spheroid_init_from_srid(gserialized_get_srid(g1), &s);
291
292 /* Set to sphere if requested */
293 if ( ! use_spheroid )
294 s.a = s.b = s.radius;
295
296 /* Return FALSE on empty arguments. */
297 if ( gserialized_is_empty(g1) || gserialized_is_empty(g2) )
298 PG_RETURN_BOOL(false);
299
300 /* Do the brute force calculation if the cached calculation doesn't tick over */
301 if (LW_FAILURE == geography_dwithin_cache(fcinfo, shared_geom1, shared_geom2, &s, tolerance, &dwithin))
302 {
303 LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
304 LWGEOM* lwgeom2 = lwgeom_from_gserialized(g2);
305 distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
306 /* Something went wrong... */
307 if ( distance < 0.0 )
308 elog(ERROR, "lwgeom_distance_spheroid returned negative!");
309 dwithin = (distance <= tolerance);
310 lwgeom_free(lwgeom1);
311 lwgeom_free(lwgeom2);
312 }
313
314 PG_RETURN_BOOL(dwithin);
315}
316
317PG_FUNCTION_INFO_V1(geography_intersects);
318Datum geography_intersects(PG_FUNCTION_ARGS)
319{
320 PG_RETURN_BOOL(CallerFInfoFunctionCall2(
321 geography_dwithin, fcinfo->flinfo, InvalidOid,
322 PG_GETARG_DATUM(0), PG_GETARG_DATUM(1)));
323}
324
325/*
326** geography_distance_tree(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
327** returns double distance in meters
328*/
329PG_FUNCTION_INFO_V1(geography_distance_tree);
330Datum geography_distance_tree(PG_FUNCTION_ARGS)
331{
332 GSERIALIZED *g1 = NULL;
333 GSERIALIZED *g2 = NULL;
334 double tolerance = 0.0;
335 double distance;
336 bool use_spheroid = true;
337 SPHEROID s;
338
339 /* Get our geometry objects loaded into memory. */
340 g1 = PG_GETARG_GSERIALIZED_P(0);
341 g2 = PG_GETARG_GSERIALIZED_P(1);
342
343 gserialized_error_if_srid_mismatch(g1, g2, __func__);
344
345 /* Return zero on empty arguments. */
346 if ( gserialized_is_empty(g1) || gserialized_is_empty(g2) )
347 {
348 PG_FREE_IF_COPY(g1, 0);
349 PG_FREE_IF_COPY(g2, 1);
350 PG_RETURN_FLOAT8(0.0);
351 }
352
353 /* Read our tolerance value. */
354 if ( PG_NARGS() > 2 && ! PG_ARGISNULL(2) )
355 tolerance = PG_GETARG_FLOAT8(2);
356
357 /* Read our calculation type. */
358 if ( PG_NARGS() > 3 && ! PG_ARGISNULL(3) )
359 use_spheroid = PG_GETARG_BOOL(3);
360
361 /* Initialize spheroid */
362 spheroid_init_from_srid(gserialized_get_srid(g1), &s);
363
364 /* Set to sphere if requested */
365 if ( ! use_spheroid )
366 s.a = s.b = s.radius;
367
368 if ( geography_tree_distance(g1, g2, &s, tolerance, &distance) == LW_FAILURE )
369 {
370 elog(ERROR, "geography_distance_tree failed!");
371 PG_RETURN_NULL();
372 }
373 /* Knock off any funny business at the nanometer level, ticket #2168 */
374 distance = round(distance * INVMINDIST) / INVMINDIST;
375
376 PG_RETURN_FLOAT8(distance);
377}
378
379
380
381/*
382** geography_dwithin_uncached(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
383** returns double distance in meters
384*/
385PG_FUNCTION_INFO_V1(geography_dwithin_uncached);
386Datum geography_dwithin_uncached(PG_FUNCTION_ARGS)
387{
388 LWGEOM *lwgeom1 = NULL;
389 LWGEOM *lwgeom2 = NULL;
390 GSERIALIZED *g1 = NULL;
391 GSERIALIZED *g2 = NULL;
392 double tolerance = 0.0;
393 double distance;
394 bool use_spheroid = true;
395 SPHEROID s;
396
397 /* Get our geometry objects loaded into memory. */
398 g1 = PG_GETARG_GSERIALIZED_P(0);
399 g2 = PG_GETARG_GSERIALIZED_P(1);
400 gserialized_error_if_srid_mismatch(g1, g2, __func__);
401
402 /* Read our tolerance value. */
403 if ( PG_NARGS() > 2 && ! PG_ARGISNULL(2) )
404 tolerance = PG_GETARG_FLOAT8(2);
405
406 /* Read our calculation type. */
407 if ( PG_NARGS() > 3 && ! PG_ARGISNULL(3) )
408 use_spheroid = PG_GETARG_BOOL(3);
409
410 /* Initialize spheroid */
411 spheroid_init_from_srid(gserialized_get_srid(g1), &s);
412
413 /* Set to sphere if requested */
414 if ( ! use_spheroid )
415 s.a = s.b = s.radius;
416
417 lwgeom1 = lwgeom_from_gserialized(g1);
418 lwgeom2 = lwgeom_from_gserialized(g2);
419
420 /* Return FALSE on empty arguments. */
421 if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
422 {
423 PG_RETURN_BOOL(false);
424 }
425
426 distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
427
428 /* Clean up */
429 lwgeom_free(lwgeom1);
430 lwgeom_free(lwgeom2);
431 PG_FREE_IF_COPY(g1, 0);
432 PG_FREE_IF_COPY(g2, 1);
433
434 /* Something went wrong... should already be eloged, return FALSE */
435 if ( distance < 0.0 )
436 {
437 elog(ERROR, "lwgeom_distance_spheroid returned negative!");
438 PG_RETURN_BOOL(false);
439 }
440
441 PG_RETURN_BOOL(distance <= tolerance);
442}
443
444
445/*
446** geography_expand(GSERIALIZED *g) returns *GSERIALIZED
447**
448** warning, this tricky little function does not expand the
449** geometry at all, just re-writes bounding box value to be
450** a bit bigger. only useful when passing the result along to
451** an index operator (&&)
452*/
453PG_FUNCTION_INFO_V1(geography_expand);
454Datum geography_expand(PG_FUNCTION_ARGS)
455{
456 GSERIALIZED *g = NULL;
457 GSERIALIZED *g_out = NULL;
458 double unit_distance, distance;
459
460 /* Get a wholly-owned pointer to the geography */
461 g = PG_GETARG_GSERIALIZED_P_COPY(0);
462
463 /* Read our distance value and normalize to unit-sphere. */
464 distance = PG_GETARG_FLOAT8(1);
465 /* Magic 1% expansion is to bridge difference between potential */
466 /* spheroidal input distance and fact that expanded box filter is */
467 /* calculated on sphere */
468 unit_distance = 1.01 * distance / WGS84_RADIUS;
469
470 /* Try the expansion */
471 g_out = gserialized_expand(g, unit_distance);
472
473 /* If the expansion fails, the return our input */
474 if ( g_out == NULL )
475 {
476 PG_RETURN_POINTER(g);
477 }
478
479 if ( g_out != g )
480 {
481 pfree(g);
482 }
483
484 PG_RETURN_POINTER(g_out);
485}
486
487/*
488** geography_area(GSERIALIZED *g)
489** returns double area in meters square
490*/
491PG_FUNCTION_INFO_V1(geography_area);
492Datum geography_area(PG_FUNCTION_ARGS)
493{
494 LWGEOM *lwgeom = NULL;
495 GSERIALIZED *g = NULL;
496 GBOX gbox;
497 double area;
498 bool use_spheroid = LW_TRUE;
499 SPHEROID s;
500
501 /* Get our geometry object loaded into memory. */
502 g = PG_GETARG_GSERIALIZED_P(0);
503
504 /* Read our calculation type */
505 use_spheroid = PG_GETARG_BOOL(1);
506
507 /* Initialize spheroid */
508 spheroid_init_from_srid(gserialized_get_srid(g), &s);
509
510 lwgeom = lwgeom_from_gserialized(g);
511
512 /* EMPTY things have no area */
513 if ( lwgeom_is_empty(lwgeom) )
514 {
515 lwgeom_free(lwgeom);
516 PG_RETURN_FLOAT8(0.0);
517 }
518
519 if ( lwgeom->bbox )
520 gbox = *(lwgeom->bbox);
521 else
522 lwgeom_calculate_gbox_geodetic(lwgeom, &gbox);
523
524#ifndef PROJ_GEODESIC
525 /* Test for cases that are currently not handled by spheroid code */
526 if ( use_spheroid )
527 {
528 /* We can't circle the poles right now */
529 if ( FP_GTEQ(gbox.zmax,1.0) || FP_LTEQ(gbox.zmin,-1.0) )
530 use_spheroid = LW_FALSE;
531 /* We can't cross the equator right now */
532 if ( gbox.zmax > 0.0 && gbox.zmin < 0.0 )
533 use_spheroid = LW_FALSE;
534 }
535#endif /* ifndef PROJ_GEODESIC */
536
537 /* User requests spherical calculation, turn our spheroid into a sphere */
538 if ( ! use_spheroid )
539 s.a = s.b = s.radius;
540
541 /* Calculate the area */
542 area = lwgeom_area_spheroid(lwgeom, &s);
543
544 /* Clean up */
545 lwgeom_free(lwgeom);
546 PG_FREE_IF_COPY(g, 0);
547
548 /* Something went wrong... */
549 if ( area < 0.0 )
550 {
551 elog(ERROR, "lwgeom_area_spher(oid) returned area < 0.0");
552 PG_RETURN_NULL();
553 }
554
555 PG_RETURN_FLOAT8(area);
556}
557
558/*
559** geography_perimeter(GSERIALIZED *g)
560** returns double perimeter in meters for area features
561*/
562PG_FUNCTION_INFO_V1(geography_perimeter);
563Datum geography_perimeter(PG_FUNCTION_ARGS)
564{
565 LWGEOM *lwgeom = NULL;
566 GSERIALIZED *g = NULL;
567 double length;
568 bool use_spheroid = LW_TRUE;
569 SPHEROID s;
570 int type;
571
572 /* Get our geometry object loaded into memory. */
573 g = PG_GETARG_GSERIALIZED_P(0);
574
575 /* Only return for area features. */
576 type = gserialized_get_type(g);
577 if ( ! (type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE) )
578 {
579 PG_RETURN_FLOAT8(0.0);
580 }
581
582 lwgeom = lwgeom_from_gserialized(g);
583
584 /* EMPTY things have no perimeter */
585 if ( lwgeom_is_empty(lwgeom) )
586 {
587 lwgeom_free(lwgeom);
588 PG_RETURN_FLOAT8(0.0);
589 }
590
591 /* Read our calculation type */
592 use_spheroid = PG_GETARG_BOOL(1);
593
594 /* Initialize spheroid */
595 spheroid_init_from_srid(gserialized_get_srid(g), &s);
596
597 /* User requests spherical calculation, turn our spheroid into a sphere */
598 if ( ! use_spheroid )
599 s.a = s.b = s.radius;
600
601 /* Calculate the length */
602 length = lwgeom_length_spheroid(lwgeom, &s);
603
604 /* Something went wrong... */
605 if ( length < 0.0 )
606 {
607 elog(ERROR, "lwgeom_length_spheroid returned length < 0.0");
608 PG_RETURN_NULL();
609 }
610
611 /* Clean up, but not all the way to the point arrays */
612 lwgeom_free(lwgeom);
613
614 PG_FREE_IF_COPY(g, 0);
615 PG_RETURN_FLOAT8(length);
616}
617
618/*
619** geography_length(GSERIALIZED *g)
620** returns double length in meters
621*/
622PG_FUNCTION_INFO_V1(geography_length);
623Datum geography_length(PG_FUNCTION_ARGS)
624{
625 LWGEOM *lwgeom = NULL;
626 GSERIALIZED *g = NULL;
627 double length;
628 bool use_spheroid = LW_TRUE;
629 SPHEROID s;
630
631 /* Get our geometry object loaded into memory. */
632 g = PG_GETARG_GSERIALIZED_P(0);
633 lwgeom = lwgeom_from_gserialized(g);
634
635 /* EMPTY things have no length */
636 if ( lwgeom_is_empty(lwgeom) || lwgeom->type == POLYGONTYPE || lwgeom->type == MULTIPOLYGONTYPE )
637 {
638 lwgeom_free(lwgeom);
639 PG_RETURN_FLOAT8(0.0);
640 }
641
642 /* Read our calculation type */
643 use_spheroid = PG_GETARG_BOOL(1);
644
645 /* Initialize spheroid */
646 spheroid_init_from_srid(gserialized_get_srid(g), &s);
647
648 /* User requests spherical calculation, turn our spheroid into a sphere */
649 if ( ! use_spheroid )
650 s.a = s.b = s.radius;
651
652 /* Calculate the length */
653 length = lwgeom_length_spheroid(lwgeom, &s);
654
655 /* Something went wrong... */
656 if ( length < 0.0 )
657 {
658 elog(ERROR, "lwgeom_length_spheroid returned length < 0.0");
659 PG_RETURN_NULL();
660 }
661
662 /* Clean up */
663 lwgeom_free(lwgeom);
664
665 PG_FREE_IF_COPY(g, 0);
666 PG_RETURN_FLOAT8(length);
667}
668
669/*
670** geography_point_outside(GSERIALIZED *g)
671** returns point outside the object
672*/
673PG_FUNCTION_INFO_V1(geography_point_outside);
674Datum geography_point_outside(PG_FUNCTION_ARGS)
675{
676 GBOX gbox;
677 LWGEOM *lwpoint = NULL;
678 POINT2D pt;
679
680 /* We need the bounding box to get an outside point for area algorithm */
681 if (gserialized_datum_get_gbox_p(PG_GETARG_DATUM(0), &gbox) == LW_FAILURE)
682 {
683 POSTGIS_DEBUG(4, "gserialized_datum_get_gbox_p returned LW_FAILURE");
684 elog(ERROR, "Error in gserialized_datum_get_gbox_p calculation.");
685 PG_RETURN_NULL();
686 }
687 POSTGIS_DEBUGF(4, "got gbox %s", gbox_to_string(&gbox));
688
689 /* Get an exterior point, based on this gbox */
690 gbox_pt_outside(&gbox, &pt);
691
692 lwpoint = (LWGEOM*) lwpoint_make2d(4326, pt.x, pt.y);
693
694 PG_RETURN_POINTER(geography_serialize(lwpoint));
695}
696
697/*
698** geography_covers(GSERIALIZED *g, GSERIALIZED *g) returns boolean
699** Only works for (multi)points and (multi)polygons currently.
700** Attempts a simple point-in-polygon test on the polygon and point.
701** Current algorithm does not distinguish between points on edge
702** and points within.
703*/
704PG_FUNCTION_INFO_V1(geography_covers);
705Datum geography_covers(PG_FUNCTION_ARGS)
706{
707 LWGEOM *lwgeom1 = NULL;
708 LWGEOM *lwgeom2 = NULL;
709 GSERIALIZED *g1 = NULL;
710 GSERIALIZED *g2 = NULL;
711 int result = LW_FALSE;
712
713 /* Get our geometry objects loaded into memory. */
714 g1 = PG_GETARG_GSERIALIZED_P(0);
715 g2 = PG_GETARG_GSERIALIZED_P(1);
716 gserialized_error_if_srid_mismatch(g1, g2, __func__);
717
718 /* Construct our working geometries */
719 lwgeom1 = lwgeom_from_gserialized(g1);
720 lwgeom2 = lwgeom_from_gserialized(g2);
721
722 /* EMPTY never intersects with another geometry */
723 if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
724 {
725 lwgeom_free(lwgeom1);
726 lwgeom_free(lwgeom2);
727 PG_FREE_IF_COPY(g1, 0);
728 PG_FREE_IF_COPY(g2, 1);
729 PG_RETURN_BOOL(false);
730 }
731
732 /* Calculate answer */
733 result = lwgeom_covers_lwgeom_sphere(lwgeom1, lwgeom2);
734
735 /* Clean up */
736 lwgeom_free(lwgeom1);
737 lwgeom_free(lwgeom2);
738 PG_FREE_IF_COPY(g1, 0);
739 PG_FREE_IF_COPY(g2, 1);
740
741 PG_RETURN_BOOL(result);
742}
743
744PG_FUNCTION_INFO_V1(geography_coveredby);
745Datum geography_coveredby(PG_FUNCTION_ARGS)
746{
747 LWGEOM *lwgeom1 = NULL;
748 LWGEOM *lwgeom2 = NULL;
749 GSERIALIZED *g1 = NULL;
750 GSERIALIZED *g2 = NULL;
751 int result = LW_FALSE;
752
753 /* Get our geometry objects loaded into memory. */
754 /* Pick them up in reverse order to covers */
755 g1 = PG_GETARG_GSERIALIZED_P(1);
756 g2 = PG_GETARG_GSERIALIZED_P(0);
757 gserialized_error_if_srid_mismatch(g1, g2, __func__);
758
759 /* Construct our working geometries */
760 lwgeom1 = lwgeom_from_gserialized(g1);
761 lwgeom2 = lwgeom_from_gserialized(g2);
762
763 /* EMPTY never intersects with another geometry */
764 if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
765 {
766 lwgeom_free(lwgeom1);
767 lwgeom_free(lwgeom2);
768 PG_FREE_IF_COPY(g1, 1);
769 PG_FREE_IF_COPY(g2, 0);
770 PG_RETURN_BOOL(false);
771 }
772
773 /* Calculate answer */
774 result = lwgeom_covers_lwgeom_sphere(lwgeom1, lwgeom2);
775
776 /* Clean up */
777 lwgeom_free(lwgeom1);
778 lwgeom_free(lwgeom2);
779 PG_FREE_IF_COPY(g1, 1);
780 PG_FREE_IF_COPY(g2, 0);
781
782 PG_RETURN_BOOL(result);
783}
784
785/*
786** geography_bestsrid(GSERIALIZED *g, GSERIALIZED *g) returns int
787** Utility function. Returns negative SRID numbers that match to the
788** numbers handled in code by the transform(lwgeom, srid) function.
789** UTM, polar stereographic and mercator as fallback. To be used
790** in wrapping existing geometry functions in SQL to provide access
791** to them in the geography module.
792*/
793PG_FUNCTION_INFO_V1(geography_bestsrid);
794Datum geography_bestsrid(PG_FUNCTION_ARGS)
795{
796 GBOX gbox, gbox1, gbox2;
797 GSERIALIZED *g1 = NULL;
798 GSERIALIZED *g2 = NULL;
799 int empty1 = LW_FALSE;
800 int empty2 = LW_FALSE;
801 double xwidth, ywidth;
802 POINT2D center;
803 LWGEOM *lwgeom;
804
805 /* Get our geometry objects loaded into memory. */
806 g1 = PG_GETARG_GSERIALIZED_P(0);
807 /* Synchronize our box types */
808 gbox1.flags = gserialized_get_lwflags(g1);
809 /* Calculate if the geometry is empty. */
810 empty1 = gserialized_is_empty(g1);
811
812 /* Convert g1 to LWGEOM type */
813 lwgeom = lwgeom_from_gserialized(g1);
814
815 /* Calculate a geocentric bounds for the objects */
816 if ( ! empty1 && gserialized_get_gbox_p(g1, &gbox1) == LW_FAILURE )
817 elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g1, &gbox1)");
818
819 POSTGIS_DEBUGF(4, "calculated gbox = %s", gbox_to_string(&gbox1));
820
821 if ( !lwgeom_isfinite(lwgeom) ) {
822 elog(ERROR, "Error in geography_bestsrid calling with infinite coordinate geographies");
823 }
824 lwgeom_free(lwgeom);
825
826 /* If we have a unique second argument, fill in all the necessary variables. */
827 if (PG_NARGS() > 1)
828 {
829 g2 = PG_GETARG_GSERIALIZED_P(1);
830 gbox2.flags = gserialized_get_lwflags(g2);
831 empty2 = gserialized_is_empty(g2);
832 if ( ! empty2 && gserialized_get_gbox_p(g2, &gbox2) == LW_FAILURE )
833 elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g2, &gbox2)");
834
835 /* Convert g2 to LWGEOM type */
836 lwgeom = lwgeom_from_gserialized(g2);
837
838 if ( !lwgeom_isfinite(lwgeom) ) {
839 elog(ERROR, "Error in geography_bestsrid calling with second arg infinite coordinate geographies");
840 }
841 lwgeom_free(lwgeom);
842 }
843 /*
844 ** If no unique second argument, copying the box for the first
845 ** argument will give us the right answer for all subsequent tests.
846 */
847 else
848 {
849 gbox = gbox2 = gbox1;
850 }
851
852 /* Both empty? We don't have an answer. */
853 if ( empty1 && empty2 )
854 PG_RETURN_NULL();
855
856 /* One empty? We can use the other argument values as infill. Otherwise merge the boxen */
857 if ( empty1 )
858 gbox = gbox2;
859 else if ( empty2 )
860 gbox = gbox1;
861 else
862 gbox_union(&gbox1, &gbox2, &gbox);
863
864 gbox_centroid(&gbox, &center);
865
866 /* Width and height in degrees */
867 xwidth = 180.0 * gbox_angular_width(&gbox) / M_PI;
868 ywidth = 180.0 * gbox_angular_height(&gbox) / M_PI;
869
870 POSTGIS_DEBUGF(2, "xwidth %g", xwidth);
871 POSTGIS_DEBUGF(2, "ywidth %g", ywidth);
872 POSTGIS_DEBUGF(2, "center POINT(%g %g)", center.x, center.y);
873
874 /* Are these data arctic? Lambert Azimuthal Equal Area North. */
875 if ( center.y > 70.0 && ywidth < 45.0 )
876 {
877 PG_RETURN_INT32(SRID_NORTH_LAMBERT);
878 }
879
880 /* Are these data antarctic? Lambert Azimuthal Equal Area South. */
881 if ( center.y < -70.0 && ywidth < 45.0 )
882 {
883 PG_RETURN_INT32(SRID_SOUTH_LAMBERT);
884 }
885
886 /*
887 ** Can we fit these data into one UTM zone?
888 ** We will assume we can push things as
889 ** far as a half zone past a zone boundary.
890 ** Note we have no handling for the date line in here.
891 */
892 if ( xwidth < 6.0 )
893 {
894 int zone = floor((center.x + 180.0) / 6.0);
895
896 if ( zone > 59 ) zone = 59;
897
898 /* Are these data below the equator? UTM South. */
899 if ( center.y < 0.0 )
900 {
901 PG_RETURN_INT32( SRID_SOUTH_UTM_START + zone );
902 }
903 /* Are these data above the equator? UTM North. */
904 else
905 {
906 PG_RETURN_INT32( SRID_NORTH_UTM_START + zone );
907 }
908 }
909
910 /*
911 ** Can we fit into a custom LAEA area? (30 degrees high, variable width)
912 ** We will allow overlap into adjoining areas, but use a slightly narrower test (25) to try
913 ** and minimize the worst case.
914 ** Again, we are hoping the dateline doesn't trip us up much
915 */
916 if ( ywidth < 25.0 )
917 {
918 int xzone = -1;
919 int yzone = 3 + floor(center.y / 30.0); /* (range of 0-5) */
920
921 /* Equatorial band, 12 zones, 30 degrees wide */
922 if ( (yzone == 2 || yzone == 3) && xwidth < 30.0 )
923 {
924 xzone = 6 + floor(center.x / 30.0);
925 }
926 /* Temperate band, 8 zones, 45 degrees wide */
927 else if ( (yzone == 1 || yzone == 4) && xwidth < 45.0 )
928 {
929 xzone = 4 + floor(center.x / 45.0);
930 }
931 /* Arctic band, 4 zones, 90 degrees wide */
932 else if ( (yzone == 0 || yzone == 5) && xwidth < 90.0 )
933 {
934 xzone = 2 + floor(center.x / 90.0);
935 }
936
937 /* Did we fit into an appropriate xzone? */
938 if ( xzone != -1 )
939 {
940 PG_RETURN_INT32(SRID_LAEA_START + 20 * yzone + xzone);
941 }
942 }
943
944 /*
945 ** Running out of options... fall-back to Mercator
946 ** and hope for the best.
947 */
948 PG_RETURN_INT32(SRID_WORLD_MERCATOR);
949
950}
951
952/*
953** geography_project(GSERIALIZED *g, distance, azimuth)
954** returns point of projection given start point,
955** azimuth in radians (bearing) and distance in meters
956*/
957PG_FUNCTION_INFO_V1(geography_project);
958Datum geography_project(PG_FUNCTION_ARGS)
959{
960 LWGEOM *lwgeom = NULL;
961 LWPOINT *lwp_projected;
962 GSERIALIZED *g = NULL;
963 GSERIALIZED *g_out = NULL;
964 double azimuth = 0.0;
965 double distance;
966 SPHEROID s;
967 uint32_t type;
968
969 /* Return NULL on NULL distance or geography */
970 if ( PG_NARGS() < 2 || PG_ARGISNULL(0) || PG_ARGISNULL(1) )
971 PG_RETURN_NULL();
972
973 /* Get our geometry object loaded into memory. */
974 g = PG_GETARG_GSERIALIZED_P(0);
975
976 /* Only return for points. */
977 type = gserialized_get_type(g);
978 if ( type != POINTTYPE )
979 {
980 elog(ERROR, "ST_Project(geography) is only valid for point inputs");
981 PG_RETURN_NULL();
982 }
983
984 distance = PG_GETARG_FLOAT8(1); /* Distance in Meters */
985 lwgeom = lwgeom_from_gserialized(g);
986
987 /* EMPTY things cannot be projected from */
988 if ( lwgeom_is_empty(lwgeom) )
989 {
990 lwgeom_free(lwgeom);
991 elog(ERROR, "ST_Project(geography) cannot project from an empty start point");
992 PG_RETURN_NULL();
993 }
994
995 if ( PG_NARGS() > 2 && ! PG_ARGISNULL(2) )
996 azimuth = PG_GETARG_FLOAT8(2); /* Azimuth in Radians */
997
998 /* Initialize spheroid */
999 spheroid_init_from_srid(gserialized_get_srid(g), &s);
1000
1001 /* Handle the zero distance case */
1002 if( FP_EQUALS(distance, 0.0) )
1003 {
1004 PG_RETURN_POINTER(g);
1005 }
1006
1007 /* Calculate the length */
1008 lwp_projected = lwgeom_project_spheroid(lwgeom_as_lwpoint(lwgeom), &s, distance, azimuth);
1009
1010 /* Something went wrong... */
1011 if ( lwp_projected == NULL )
1012 {
1013 elog(ERROR, "lwgeom_project_spheroid returned null");
1014 PG_RETURN_NULL();
1015 }
1016
1017 /* Clean up, but not all the way to the point arrays */
1018 lwgeom_free(lwgeom);
1019 g_out = geography_serialize(lwpoint_as_lwgeom(lwp_projected));
1020 lwpoint_free(lwp_projected);
1021
1022 PG_FREE_IF_COPY(g, 0);
1023 PG_RETURN_POINTER(g_out);
1024}
1025
1026/*
1027** geography_project_geography(geog1, geog2, distance, azimuth)
1028** returns point of projection given from/pt points,
1029** azimuth in radians (bearing) and distance in meters
1030*/
1031PG_FUNCTION_INFO_V1(geography_project_geography);
1032Datum geography_project_geography(PG_FUNCTION_ARGS)
1033{
1034 LWGEOM *lwgeom1, *lwgeom2;
1035 LWPOINT *lwp1, *lwp2, *lwp3;
1036 GSERIALIZED *g1, *g2, *g3;
1037 double distance;
1038 SPHEROID s;
1039
1040 /* Get our geometry object loaded into memory. */
1041 g1 = PG_GETARG_GSERIALIZED_P(0);
1042 g2 = PG_GETARG_GSERIALIZED_P(1);
1043
1044 if ( gserialized_get_type(g1) != POINTTYPE ||
1045 gserialized_get_type(g2) != POINTTYPE )
1046 {
1047 elog(ERROR, "ST_Project(geography) is only valid for point inputs");
1048 PG_RETURN_NULL();
1049 }
1050
1051 distance = PG_GETARG_FLOAT8(2); /* Distance in meters */
1052
1053 /* Handle the zero distance case */
1054 if ( FP_EQUALS(distance, 0.0) )
1055 {
1056 PG_RETURN_POINTER(g2);
1057 }
1058
1059 lwgeom1 = lwgeom_from_gserialized(g1);
1060 lwgeom2 = lwgeom_from_gserialized(g2);
1061
1062 /* EMPTY things cannot be projected from */
1063 if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
1064 {
1065 lwgeom_free(lwgeom1);
1066 lwgeom_free(lwgeom2);
1067 elog(ERROR, "ST_Project(geography) cannot project from an empty point");
1068 PG_RETURN_NULL();
1069 }
1070
1071 /* Initialize spheroid */
1072 spheroid_init_from_srid(lwgeom_get_srid(lwgeom1), &s);
1073
1074 /* Calculate the length */
1075 lwp1 = lwgeom_as_lwpoint(lwgeom1);
1076 lwp2 = lwgeom_as_lwpoint(lwgeom2);
1077 lwp3 = lwgeom_project_spheroid_lwpoint(lwp1, lwp2, &s, distance);
1078
1079 /* Something went wrong... */
1080 if ( lwp3 == NULL )
1081 {
1082 elog(ERROR, "lwgeom_project_spheroid_lwpoint returned null");
1083 PG_RETURN_NULL();
1084 }
1085
1086 /* Clean up, but not all the way to the point arrays */
1087 lwgeom_free(lwgeom1);
1088 lwgeom_free(lwgeom2);
1089 g3 = geography_serialize(lwpoint_as_lwgeom(lwp3));
1090 lwpoint_free(lwp3);
1091
1092 PG_FREE_IF_COPY(g1, 0);
1093 PG_FREE_IF_COPY(g2, 1);
1094 PG_RETURN_POINTER(g3);
1095}
1096
1097
1098/*
1099** geography_azimuth(GSERIALIZED *g1, GSERIALIZED *g2)
1100** returns direction between points (north = 0)
1101** azimuth (bearing) and distance
1102*/
1103PG_FUNCTION_INFO_V1(geography_azimuth);
1104Datum geography_azimuth(PG_FUNCTION_ARGS)
1105{
1106 GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
1107 GSERIALIZED *g2 = PG_GETARG_GSERIALIZED_P(1);
1108 LWGEOM *lwgeom1 = NULL;
1109 LWGEOM *lwgeom2 = NULL;
1110 double azimuth;
1111 SPHEROID s;
1112 uint32_t type1, type2;
1113
1114 /* Only return for points. */
1115 type1 = gserialized_get_type(g1);
1116 type2 = gserialized_get_type(g2);
1117 if ( type1 != POINTTYPE || type2 != POINTTYPE )
1118 {
1119 elog(ERROR, "ST_Azimuth(geography, geography) is only valid for point inputs");
1120 PG_RETURN_NULL();
1121 }
1122
1123 lwgeom1 = lwgeom_from_gserialized(g1);
1124 lwgeom2 = lwgeom_from_gserialized(g2);
1125
1126 /* EMPTY things cannot be used */
1127 if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
1128 {
1129 lwgeom_free(lwgeom1);
1130 lwgeom_free(lwgeom2);
1131 elog(ERROR, "ST_Azimuth(geography, geography) cannot work with empty points");
1132 PG_RETURN_NULL();
1133 }
1134
1135 /* Initialize spheroid */
1136 spheroid_init_from_srid(gserialized_get_srid(g1), &s);
1137
1138 /* Calculate the direction */
1139 azimuth = lwgeom_azumith_spheroid(lwgeom_as_lwpoint(lwgeom1), lwgeom_as_lwpoint(lwgeom2), &s);
1140
1141 /* Clean up */
1142 lwgeom_free(lwgeom1);
1143 lwgeom_free(lwgeom2);
1144
1145 PG_FREE_IF_COPY(g1, 0);
1146 PG_FREE_IF_COPY(g2, 1);
1147
1148 /* Return NULL for unknown (same point) azimuth */
1149 if( !isfinite(azimuth) )
1150 {
1151 PG_RETURN_NULL();
1152 }
1153
1154 PG_RETURN_FLOAT8(azimuth);
1155}
1156
1157
1158
1163PG_FUNCTION_INFO_V1(geography_segmentize);
1164Datum geography_segmentize(PG_FUNCTION_ARGS)
1165{
1166 LWGEOM *lwgeom1 = NULL;
1167 LWGEOM *lwgeom2 = NULL;
1168 GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
1169 GSERIALIZED *g2 = NULL;
1170 double max_seg_length = PG_GETARG_FLOAT8(1) / WGS84_RADIUS;
1171 uint32_t type1 = gserialized_get_type(g1);
1172
1173 /* We can't densify points or points, reflect them back */
1174 if ( type1 == POINTTYPE || type1 == MULTIPOINTTYPE || gserialized_is_empty(g1) )
1175 PG_RETURN_POINTER(g1);
1176
1177 /* Deserialize */
1178 lwgeom1 = lwgeom_from_gserialized(g1);
1179
1180 /* Calculate the densified geometry */
1181 lwgeom2 = lwgeom_segmentize_sphere(lwgeom1, max_seg_length);
1182
1183 /*
1184 ** Set the geodetic flag so subsequent
1185 ** functions do the right thing.
1186 */
1187 lwgeom_set_geodetic(lwgeom2, true);
1188
1189 /* Recalculate the boxes after re-setting the geodetic bit */
1190 lwgeom_drop_bbox(lwgeom2);
1191
1192 /* We are trusting geography_serialize will add a box if needed */
1193 g2 = geography_serialize(lwgeom2);
1194
1195 /* Clean up */
1196 lwgeom_free(lwgeom1);
1197 lwgeom_free(lwgeom2);
1198 PG_FREE_IF_COPY(g1, 0);
1199
1200 PG_RETURN_POINTER(g2);
1201}
1202
1203
1204/********************************************************************************/
1205
1210PG_FUNCTION_INFO_V1(geography_line_substring);
1211Datum geography_line_substring(PG_FUNCTION_ARGS)
1212{
1213 GSERIALIZED *gs = PG_GETARG_GSERIALIZED_P(0);
1214 double from_fraction = PG_GETARG_FLOAT8(1);
1215 double to_fraction = PG_GETARG_FLOAT8(2);
1216 LWLINE *lwline;
1217 LWGEOM *lwresult;
1218 SPHEROID s;
1219 GSERIALIZED *result;
1220 bool use_spheroid = true;
1221
1222 if ( PG_NARGS() > 3 && ! PG_ARGISNULL(3) )
1223 use_spheroid = PG_GETARG_BOOL(3);
1224
1225 /* Return NULL on empty argument. */
1226 if ( gserialized_is_empty(gs) )
1227 {
1228 PG_FREE_IF_COPY(gs, 0);
1229 PG_RETURN_NULL();
1230 }
1231
1232 if ( from_fraction < 0 || from_fraction > 1 )
1233 {
1234 elog(ERROR,"%s: second argument is not within [0,1]", __func__);
1235 PG_FREE_IF_COPY(gs, 0);
1236 PG_RETURN_NULL();
1237 }
1238 if ( to_fraction < 0 || to_fraction > 1 )
1239 {
1240 elog(ERROR,"%s: argument arg is not within [0,1]", __func__);
1241 PG_FREE_IF_COPY(gs, 0);
1242 PG_RETURN_NULL();
1243 }
1244 if ( from_fraction > to_fraction )
1245 {
1246 elog(ERROR, "%s: second argument must be smaller than third argument", __func__);
1247 PG_RETURN_NULL();
1248 }
1249
1250 lwline = lwgeom_as_lwline(lwgeom_from_gserialized(gs));
1251 if ( !lwline )
1252 {
1253 elog(ERROR,"%s: first argument is not a line", __func__);
1254 PG_FREE_IF_COPY(gs, 0);
1255 PG_RETURN_NULL();
1256 }
1257
1258 /* Initialize spheroid */
1259 spheroid_init_from_srid(gserialized_get_srid(gs), &s);
1260 /* Set to sphere if requested */
1261 if ( ! use_spheroid )
1262 s.a = s.b = s.radius;
1263
1264 lwresult = geography_substring(lwline, &s,
1265 from_fraction, to_fraction, FP_TOLERANCE);
1266
1267 lwline_free(lwline);
1268 PG_FREE_IF_COPY(gs, 0);
1269 lwgeom_set_geodetic(lwresult, true);
1270 result = geography_serialize(lwresult);
1271 lwgeom_free(lwresult);
1272
1273 PG_RETURN_POINTER(result);
1274}
1275
1276
1285PG_FUNCTION_INFO_V1(geography_line_interpolate_point);
1286Datum geography_line_interpolate_point(PG_FUNCTION_ARGS)
1287{
1288 GSERIALIZED *gs = PG_GETARG_GSERIALIZED_P(0);
1289 double distance_fraction = PG_GETARG_FLOAT8(1);
1290 /* Read calculation type */
1291 bool use_spheroid = PG_GETARG_BOOL(2);
1292 /* Read repeat mode */
1293 bool repeat = (PG_NARGS() > 3) && PG_GETARG_BOOL(3);
1294 LWLINE* lwline;
1295 LWGEOM* lwresult;
1296 SPHEROID s;
1297 GSERIALIZED *result;
1298
1299 /* Return NULL on empty argument. */
1300 if ( gserialized_is_empty(gs) )
1301 {
1302 PG_FREE_IF_COPY(gs, 0);
1303 PG_RETURN_NULL();
1304 }
1305
1306 if ( distance_fraction < 0 || distance_fraction > 1 )
1307 {
1308 elog(ERROR,"%s: second arg is not within [0,1]", __func__);
1309 PG_FREE_IF_COPY(gs, 0);
1310 PG_RETURN_NULL();
1311 }
1312
1313 lwline = lwgeom_as_lwline(lwgeom_from_gserialized(gs));
1314 if ( !lwline )
1315 {
1316 elog(ERROR,"%s: first arg is not a line", __func__);
1317 PG_FREE_IF_COPY(gs, 0);
1318 PG_RETURN_NULL();
1319 }
1320
1321 /* Initialize spheroid */
1322 spheroid_init_from_srid(gserialized_get_srid(gs), &s);
1323
1324 /* Set to sphere if requested */
1325 if ( ! use_spheroid )
1326 s.a = s.b = s.radius;
1327
1328 lwresult = geography_interpolate_points(lwline, distance_fraction, &s, repeat);
1329
1330 lwgeom_free(lwline_as_lwgeom(lwline));
1331 PG_FREE_IF_COPY(gs, 0);
1332
1333 lwgeom_set_geodetic(lwresult, true);
1334 result = geography_serialize(lwresult);
1335 lwgeom_free(lwresult);
1336
1337 PG_RETURN_POINTER(result);
1338}
1339
1340
1345PG_FUNCTION_INFO_V1(geography_line_locate_point);
1346Datum geography_line_locate_point(PG_FUNCTION_ARGS)
1347{
1348 GSERIALIZED *gs1 = PG_GETARG_GSERIALIZED_P(0);
1349 GSERIALIZED *gs2 = PG_GETARG_GSERIALIZED_P(1);
1350 bool use_spheroid = PG_GETARG_BOOL(2);
1351 double tolerance = FP_TOLERANCE;
1352 SPHEROID s;
1353 LWLINE *lwline;
1354 LWPOINT *lwpoint;
1355 POINTARRAY *pa;
1356 POINT4D p, p_proj;
1357 double ret;
1358
1359 gserialized_error_if_srid_mismatch(gs1, gs2, __func__);
1360
1361 /* Return NULL on empty argument. */
1362 if ( gserialized_is_empty(gs1) || gserialized_is_empty(gs2))
1363 {
1364 PG_FREE_IF_COPY(gs1, 0);
1365 PG_FREE_IF_COPY(gs2, 1);
1366 PG_RETURN_NULL();
1367 }
1368
1369 if ( gserialized_get_type(gs1) != LINETYPE )
1370 {
1371 elog(ERROR,"%s: 1st arg is not a line", __func__);
1372 PG_RETURN_NULL();
1373 }
1374 if ( gserialized_get_type(gs2) != POINTTYPE )
1375 {
1376 elog(ERROR,"%s: 2nd arg is not a point", __func__);
1377 PG_RETURN_NULL();
1378 }
1379
1380 /* Set to sphere if requested */
1381 if ( ! use_spheroid ) {
1382 s.a = s.b = s.radius;
1383 }
1384 else {
1385 /* Initialize spheroid */
1386 spheroid_init_from_srid(gserialized_get_srid(gs1), &s);
1387 }
1388
1389 lwline = lwgeom_as_lwline(lwgeom_from_gserialized(gs1));
1390 lwpoint = lwgeom_as_lwpoint(lwgeom_from_gserialized(gs2));
1391
1392 pa = lwline->points;
1393 lwpoint_getPoint4d_p(lwpoint, &p);
1394
1395 ret = ptarray_locate_point_spheroid(pa, &p, &s, tolerance, NULL, &p_proj);
1396
1397 PG_RETURN_FLOAT8(ret);
1398}
1399
1400
1406PG_FUNCTION_INFO_V1(geography_closestpoint);
1407Datum geography_closestpoint(PG_FUNCTION_ARGS)
1408{
1409 GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
1410 GSERIALIZED *g2 = PG_GETARG_GSERIALIZED_P(1);
1411 LWGEOM *point, *lwg1, *lwg2;
1412 GSERIALIZED *result;
1413
1414 gserialized_error_if_srid_mismatch(g1, g2, __func__);
1415
1416 lwg1 = lwgeom_from_gserialized(g1);
1417 lwg2 = lwgeom_from_gserialized(g2);
1418
1419 /* Return NULL on empty/bad arguments. */
1420 if ( !lwg1 || !lwg2 || lwgeom_is_empty(lwg1) || lwgeom_is_empty(lwg2) )
1421 {
1422 PG_FREE_IF_COPY(g1, 0);
1423 PG_FREE_IF_COPY(g2, 1);
1424 PG_RETURN_NULL();
1425 }
1426
1427 point = geography_tree_closestpoint(lwg1, lwg2, FP_TOLERANCE);
1428 result = geography_serialize(point);
1429 lwgeom_free(point);
1430 lwgeom_free(lwg1);
1431 lwgeom_free(lwg2);
1432
1433 PG_FREE_IF_COPY(g1, 0);
1434 PG_FREE_IF_COPY(g2, 1);
1435 PG_RETURN_POINTER(result);
1436}
1437
1442PG_FUNCTION_INFO_V1(geography_shortestline);
1443Datum geography_shortestline(PG_FUNCTION_ARGS)
1444{
1445 GSERIALIZED* g1 = PG_GETARG_GSERIALIZED_P(0);
1446 GSERIALIZED* g2 = PG_GETARG_GSERIALIZED_P(1);
1447 bool use_spheroid = PG_GETARG_BOOL(2);
1448 LWGEOM *line, *lwgeom1, *lwgeom2;
1449 GSERIALIZED* result;
1450 SPHEROID s;
1451
1452 gserialized_error_if_srid_mismatch(g1, g2, __func__);
1453
1454 lwgeom1 = lwgeom_from_gserialized(g1);
1455 lwgeom2 = lwgeom_from_gserialized(g2);
1456
1457 /* Return NULL on empty/bad arguments. */
1458 if ( !lwgeom1 || !lwgeom2 || lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
1459 {
1460 PG_FREE_IF_COPY(g1, 0);
1461 PG_FREE_IF_COPY(g2, 1);
1462 PG_RETURN_NULL();
1463 }
1464
1465 /* Initialize spheroid */
1466 spheroid_init_from_srid(gserialized_get_srid(g1), &s);
1467
1468 /* Set to sphere if requested */
1469 if ( ! use_spheroid )
1470 s.a = s.b = s.radius;
1471
1472 line = geography_tree_shortestline(lwgeom1, lwgeom2, FP_TOLERANCE, &s);
1473
1474 if (lwgeom_is_empty(line))
1475 PG_RETURN_NULL();
1476
1477 result = geography_serialize(line);
1478 lwgeom_free(line);
1479
1480 PG_FREE_IF_COPY(g1, 0);
1481 PG_FREE_IF_COPY(g2, 1);
1482 PG_RETURN_POINTER(result);
1483}
s
char * s
Definition cu_in_wkt.c:23
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition cu_print.c:267
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:404
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:623
geography_area
Datum geography_area(PG_FUNCTION_ARGS)
Definition geography_measurement.c:492
geography_expand
Datum geography_expand(PG_FUNCTION_ARGS)
Definition geography_measurement.c:454
PG_FUNCTION_INFO_V1
PG_FUNCTION_INFO_V1(geography_distance_knn)
geography_azimuth
Datum geography_azimuth(PG_FUNCTION_ARGS)
Definition geography_measurement.c:1104
geography_line_substring
Datum geography_line_substring(PG_FUNCTION_ARGS)
Definition geography_measurement.c:1211
geography_line_locate_point
Datum geography_line_locate_point(PG_FUNCTION_ARGS)
Definition geography_measurement.c:1346
geography_shortestline
Datum geography_shortestline(PG_FUNCTION_ARGS)
Definition geography_measurement.c:1443
geography_closestpoint
Datum geography_closestpoint(PG_FUNCTION_ARGS)
Definition geography_measurement.c:1407
geography_distance_tree
Datum geography_distance_tree(PG_FUNCTION_ARGS)
Definition geography_measurement.c:330
geography_project_geography
Datum geography_project_geography(PG_FUNCTION_ARGS)
Definition geography_measurement.c:1032
geography_point_outside
Datum geography_point_outside(PG_FUNCTION_ARGS)
Definition geography_measurement.c:674
geography_distance_uncached
Datum geography_distance_uncached(PG_FUNCTION_ARGS)
Definition geography_measurement.c:140
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:1164
geography_covers
Datum geography_covers(PG_FUNCTION_ARGS)
Definition geography_measurement.c:705
geography_intersects
Datum geography_intersects(PG_FUNCTION_ARGS)
Definition geography_measurement.c:318
geography_dwithin
Datum geography_dwithin(PG_FUNCTION_ARGS)
Definition geography_measurement.c:267
geography_coveredby
Datum geography_coveredby(PG_FUNCTION_ARGS)
Definition geography_measurement.c:745
geography_perimeter
Datum geography_perimeter(PG_FUNCTION_ARGS)
Definition geography_measurement.c:563
geography_bestsrid
Datum geography_bestsrid(PG_FUNCTION_ARGS)
Definition geography_measurement.c:794
geography_dwithin_uncached
Datum geography_dwithin_uncached(PG_FUNCTION_ARGS)
Definition geography_measurement.c:386
geography_project
Datum geography_project(PG_FUNCTION_ARGS)
Definition geography_measurement.c:958
INVMINDIST
#define INVMINDIST
Definition geography_measurement.c:48
geography_line_interpolate_point
Datum geography_line_interpolate_point(PG_FUNCTION_ARGS)
Definition geography_measurement.c:1286
geography_distance_knn
Datum geography_distance_knn(PG_FUNCTION_ARGS)
Definition geography_measurement.c:77
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:287
geography_dwithin_cache
int geography_dwithin_cache(FunctionCallInfo fcinfo, SHARED_GSERIALIZED *g1, SHARED_GSERIALIZED *g2, const SPHEROID *s, double tolerance, int *dwithin)
Definition geography_measurement_trees.c:264
geography_distance_cache
int geography_distance_cache(FunctionCallInfo fcinfo, SHARED_GSERIALIZED *g1, SHARED_GSERIALIZED *g2, const SPHEROID *s, double *distance)
Definition geography_measurement_trees.c:254
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:432
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:155
lwgeom_from_gserialized
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Definition gserialized.c:268
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:94
gserialized_is_empty
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition gserialized.c:181
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:118
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:47
gserialized_datum_get_gbox_p
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...
Definition gserialized_gist_2d.c:514
lwpoint_getPoint4d_p
int lwpoint_getPoint4d_p(const LWPOINT *point, POINT4D *out)
Definition lwpoint.c:57
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:992
lwpoint_as_lwgeom
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition lwgeom.c:372
LW_FALSE
#define LW_FALSE
Definition liblwgeom.h:94
COLLECTIONTYPE
#define COLLECTIONTYPE
Definition liblwgeom.h:108
lwgeom_get_srid
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition lwgeom.c:955
ptarray_locate_point_spheroid
double ptarray_locate_point_spheroid(const POINTARRAY *pa, const POINT4D *p4d, const SPHEROID *s, double tolerance, double *mindistout, POINT4D *proj4d)
Locate a point along the point array defining a geographic line.
Definition lwgeodetic_measures.c:362
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:2066
LW_FAILURE
#define LW_FAILURE
Definition liblwgeom.h:96
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition lwgeom.c:1246
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:1506
lwgeom_project_spheroid_lwpoint
LWPOINT * lwgeom_project_spheroid_lwpoint(const LWPOINT *from, const LWPOINT *to, const SPHEROID *spheroid, double distance)
Calculate the location of a point on a spheroid, give a start point, end point and distance.
Definition lwgeodetic.c:2016
LINETYPE
#define LINETYPE
Definition liblwgeom.h:103
geography_substring
LWGEOM * geography_substring(const LWLINE *line, const SPHEROID *s, double from, double to, double tolerance)
Return the part of a line between two fractional locations.
Definition lwgeodetic_measures.c:104
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:2876
WGS84_RADIUS
#define WGS84_RADIUS
Definition liblwgeom.h:148
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:105
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:2292
POINTTYPE
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition liblwgeom.h:102
lwgeom_drop_bbox
void lwgeom_drop_bbox(LWGEOM *lwgeom)
Call this function to drop BBOX and SRID from LWGEOM.
Definition lwgeom.c:710
lwgeom_isfinite
int lwgeom_isfinite(const LWGEOM *lwgeom)
Check if a LWGEOM has any non-finite (NaN or Inf) coordinates.
Definition lwgeom.c:2835
MULTIPOLYGONTYPE
#define MULTIPOLYGONTYPE
Definition liblwgeom.h:107
lwpoint_make2d
LWPOINT * lwpoint_make2d(int32_t srid, double x, double y)
Definition lwpoint.c:163
POLYGONTYPE
#define POLYGONTYPE
Definition liblwgeom.h:104
lwline_as_lwgeom
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition lwgeom.c:367
lwgeom_as_lwline
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition lwgeom.c:207
geography_interpolate_points
LWGEOM * geography_interpolate_points(const LWLINE *line, double length_fraction, const SPHEROID *s, char repeat)
Interpolate a point along a geographic line.
Definition lwgeodetic_measures.c:257
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:742
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:3145
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:1967
LW_TRUE
#define LW_TRUE
Return types for functions with status returns.
Definition liblwgeom.h:93
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:2032
lwline_free
void lwline_free(LWLINE *line)
Definition lwline.c:67
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:1697
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
geography_tree_closestpoint
LWGEOM * geography_tree_closestpoint(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, double threshold)
Definition lwgeodetic_tree.c:1035
geography_tree_shortestline
LWGEOM * geography_tree_shortestline(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, double threshold, const SPHEROID *spheroid)
Definition lwgeodetic_tree.c:1069
lwgeom_as_lwpoint
static LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition lwinline.h:127
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:199
distance
static double distance(double x1, double y1, double x2, double y2)
Definition lwtree.c:1032
GBOX::zmax
double zmax
Definition liblwgeom.h:359
GBOX::zmin
double zmin
Definition liblwgeom.h:358
GBOX::flags
lwflags_t flags
Definition liblwgeom.h:353
LWGEOM::type
uint8_t type
Definition liblwgeom.h:462
LWGEOM::bbox
GBOX * bbox
Definition liblwgeom.h:458
LWLINE::points
POINTARRAY * points
Definition liblwgeom.h:483
POINT2D::y
double y
Definition liblwgeom.h:390
POINT2D::x
double x
Definition liblwgeom.h:390