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