PostGIS  2.3.8dev-r@@SVN_REVISION@@
lwgeodetic.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 2009 Paul Ramsey <pramsey@cleverelephant.ca>
22  * Copyright 2009 David Skea <David.Skea@gov.bc.ca>
23  *
24  **********************************************************************/
25 
26 
27 #include "liblwgeom_internal.h"
28 #include "lwgeodetic.h"
29 #include "lwgeom_log.h"
30 
37 
41 double longitude_radians_normalize(double lon)
42 {
43  if ( lon == -1.0 * M_PI )
44  return M_PI;
45  if ( lon == -2.0 * M_PI )
46  return 0.0;
47 
48  if ( lon > 2.0 * M_PI )
49  lon = remainder(lon, 2.0 * M_PI);
50 
51  if ( lon < -2.0 * M_PI )
52  lon = remainder(lon, -2.0 * M_PI);
53 
54  if ( lon > M_PI )
55  lon = -2.0 * M_PI + lon;
56 
57  if ( lon < -1.0 * M_PI )
58  lon = 2.0 * M_PI + lon;
59 
60  if ( lon == -2.0 * M_PI )
61  lon *= -1.0;
62 
63  return lon;
64 }
65 
69 double latitude_radians_normalize(double lat)
70 {
71 
72  if ( lat > 2.0 * M_PI )
73  lat = remainder(lat, 2.0 * M_PI);
74 
75  if ( lat < -2.0 * M_PI )
76  lat = remainder(lat, -2.0 * M_PI);
77 
78  if ( lat > M_PI )
79  lat = M_PI - lat;
80 
81  if ( lat < -1.0 * M_PI )
82  lat = -1.0 * M_PI - lat;
83 
84  if ( lat > M_PI_2 )
85  lat = M_PI - lat;
86 
87  if ( lat < -1.0 * M_PI_2 )
88  lat = -1.0 * M_PI - lat;
89 
90  return lat;
91 }
92 
97 double longitude_degrees_normalize(double lon)
98 {
99  if ( lon > 360.0 )
100  lon = remainder(lon, 360.0);
101 
102  if ( lon < -360.0 )
103  lon = remainder(lon, -360.0);
104 
105  if ( lon > 180.0 )
106  lon = -360.0 + lon;
107 
108  if ( lon < -180.0 )
109  lon = 360 + lon;
110 
111  if ( lon == -180.0 )
112  return 180.0;
113 
114  if ( lon == -360.0 )
115  return 0.0;
116 
117  return lon;
118 }
119 
124 double latitude_degrees_normalize(double lat)
125 {
126 
127  if ( lat > 360.0 )
128  lat = remainder(lat, 360.0);
129 
130  if ( lat < -360.0 )
131  lat = remainder(lat, -360.0);
132 
133  if ( lat > 180.0 )
134  lat = 180.0 - lat;
135 
136  if ( lat < -180.0 )
137  lat = -180.0 - lat;
138 
139  if ( lat > 90.0 )
140  lat = 180.0 - lat;
141 
142  if ( lat < -90.0 )
143  lat = -180.0 - lat;
144 
145  return lat;
146 }
147 
151 void point_shift(GEOGRAPHIC_POINT *p, double shift)
152 {
153  double lon = p->lon + shift;
154  if ( lon > M_PI )
155  p->lon = -1.0 * M_PI + (lon - M_PI);
156  else
157  p->lon = lon;
158  return;
159 }
160 
162 {
163  return FP_EQUALS(g1->lat, g2->lat) && FP_EQUALS(g1->lon, g2->lon);
164 }
165 
171 void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
172 {
175 }
176 
178 double
180 {
181  double d[6];
182  int i;
183  double zmin = FLT_MAX;
184  double zmax = -1 * FLT_MAX;
185  POINT3D pt;
186 
187  /* Take a copy of the box corners so we can treat them as a list */
188  /* Elements are xmin, xmax, ymin, ymax, zmin, zmax */
189  memcpy(d, &(gbox->xmin), 6*sizeof(double));
190 
191  /* Generate all 8 corner vectors of the box */
192  for ( i = 0; i < 8; i++ )
193  {
194  pt.x = d[i / 4];
195  pt.y = d[2 + (i % 4) / 2];
196  pt.z = d[4 + (i % 2)];
197  normalize(&pt);
198  if ( pt.z < zmin ) zmin = pt.z;
199  if ( pt.z > zmax ) zmax = pt.z;
200  }
201  return asin(zmax) - asin(zmin);
202 }
203 
205 double
207 {
208  double d[6];
209  int i, j;
210  POINT3D pt[3];
211  double maxangle;
212  double magnitude;
213 
214  /* Take a copy of the box corners so we can treat them as a list */
215  /* Elements are xmin, xmax, ymin, ymax, zmin, zmax */
216  memcpy(d, &(gbox->xmin), 6*sizeof(double));
217 
218  /* Start with the bottom corner */
219  pt[0].x = gbox->xmin;
220  pt[0].y = gbox->ymin;
221  magnitude = sqrt(pt[0].x*pt[0].x + pt[0].y*pt[0].y);
222  pt[0].x /= magnitude;
223  pt[0].y /= magnitude;
224 
225  /* Generate all 8 corner vectors of the box */
226  /* Find the vector furthest from our seed vector */
227  for ( j = 0; j < 2; j++ )
228  {
229  maxangle = -1 * FLT_MAX;
230  for ( i = 0; i < 4; i++ )
231  {
232  double angle, dotprod;
233  POINT3D pt_n;
234 
235  pt_n.x = d[i / 2];
236  pt_n.y = d[2 + (i % 2)];
237  magnitude = sqrt(pt_n.x*pt_n.x + pt_n.y*pt_n.y);
238  pt_n.x /= magnitude;
239  pt_n.y /= magnitude;
240  pt_n.z = 0.0;
241 
242  dotprod = pt_n.x*pt[j].x + pt_n.y*pt[j].y;
243  angle = acos(dotprod > 1.0 ? 1.0 : dotprod);
244  if ( angle > maxangle )
245  {
246  pt[j+1] = pt_n;
247  maxangle = angle;
248  }
249  }
250  }
251 
252  /* Return the distance between the two furthest vectors */
253  return maxangle;
254 }
255 
257 int
258 gbox_centroid(const GBOX* gbox, POINT2D* out)
259 {
260  double d[6];
262  POINT3D pt;
263  int i;
264 
265  /* Take a copy of the box corners so we can treat them as a list */
266  /* Elements are xmin, xmax, ymin, ymax, zmin, zmax */
267  memcpy(d, &(gbox->xmin), 6*sizeof(double));
268 
269  /* Zero out our return vector */
270  pt.x = pt.y = pt.z = 0.0;
271 
272  for ( i = 0; i < 8; i++ )
273  {
274  POINT3D pt_n;
275 
276  pt_n.x = d[i / 4];
277  pt_n.y = d[2 + ((i % 4) / 2)];
278  pt_n.z = d[4 + (i % 2)];
279  normalize(&pt_n);
280 
281  pt.x += pt_n.x;
282  pt.y += pt_n.y;
283  pt.z += pt_n.z;
284  }
285 
286  pt.x /= 8.0;
287  pt.y /= 8.0;
288  pt.z /= 8.0;
289  normalize(&pt);
290 
291  cart2geog(&pt, &g);
294 
295  return LW_SUCCESS;
296 }
297 
307 static int gbox_check_poles(GBOX *gbox)
308 {
309  int rv = LW_FALSE;
310  LWDEBUG(4, "checking poles");
311  LWDEBUGF(4, "gbox %s", gbox_to_string(gbox));
312  /* Z axis */
313  if ( gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
314  gbox->ymin < 0.0 && gbox->ymax > 0.0 )
315  {
316  if ( (gbox->zmin + gbox->zmax) > 0.0 )
317  {
318  LWDEBUG(4, "enclosed positive z axis");
319  gbox->zmax = 1.0;
320  }
321  else
322  {
323  LWDEBUG(4, "enclosed negative z axis");
324  gbox->zmin = -1.0;
325  }
326  rv = LW_TRUE;
327  }
328 
329  /* Y axis */
330  if ( gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
331  gbox->zmin < 0.0 && gbox->zmax > 0.0 )
332  {
333  if ( gbox->ymin + gbox->ymax > 0.0 )
334  {
335  LWDEBUG(4, "enclosed positive y axis");
336  gbox->ymax = 1.0;
337  }
338  else
339  {
340  LWDEBUG(4, "enclosed negative y axis");
341  gbox->ymin = -1.0;
342  }
343  rv = LW_TRUE;
344  }
345 
346  /* X axis */
347  if ( gbox->ymin < 0.0 && gbox->ymax > 0.0 &&
348  gbox->zmin < 0.0 && gbox->zmax > 0.0 )
349  {
350  if ( gbox->xmin + gbox->xmax > 0.0 )
351  {
352  LWDEBUG(4, "enclosed positive x axis");
353  gbox->xmax = 1.0;
354  }
355  else
356  {
357  LWDEBUG(4, "enclosed negative x axis");
358  gbox->xmin = -1.0;
359  }
360  rv = LW_TRUE;
361  }
362 
363  return rv;
364 }
365 
370 {
371  p->x = cos(g->lat) * cos(g->lon);
372  p->y = cos(g->lat) * sin(g->lon);
373  p->z = sin(g->lat);
374 }
375 
380 {
381  g->lon = atan2(p->y, p->x);
382  g->lat = asin(p->z);
383 }
384 
388 void ll2cart(const POINT2D *g, POINT3D *p)
389 {
390  double x_rad = M_PI * g->x / 180.0;
391  double y_rad = M_PI * g->y / 180.0;
392  double cos_y_rad = cos(y_rad);
393  p->x = cos_y_rad * cos(x_rad);
394  p->y = cos_y_rad * sin(x_rad);
395  p->z = sin(y_rad);
396 }
397 
411 static double dot_product(const POINT3D *p1, const POINT3D *p2)
412 {
413  return (p1->x*p2->x) + (p1->y*p2->y) + (p1->z*p2->z);
414 }
415 
419 static void cross_product(const POINT3D *a, const POINT3D *b, POINT3D *n)
420 {
421  n->x = a->y * b->z - a->z * b->y;
422  n->y = a->z * b->x - a->x * b->z;
423  n->z = a->x * b->y - a->y * b->x;
424  return;
425 }
426 
430 void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n)
431 {
432  n->x = a->x + b->x;
433  n->y = a->y + b->y;
434  n->z = a->z + b->z;
435  return;
436 }
437 
441 static void vector_difference(const POINT3D *a, const POINT3D *b, POINT3D *n)
442 {
443  n->x = a->x - b->x;
444  n->y = a->y - b->y;
445  n->z = a->z - b->z;
446  return;
447 }
448 
452 static void vector_scale(POINT3D *n, double scale)
453 {
454  n->x *= scale;
455  n->y *= scale;
456  n->z *= scale;
457  return;
458 }
459 
460 /*
461 * static inline double vector_magnitude(const POINT3D* v)
462 * {
463 * return sqrt(v->x*v->x + v->y*v->y + v->z*v->z);
464 * }
465 */
466 
470 double vector_angle(const POINT3D* v1, const POINT3D* v2)
471 {
472  POINT3D v3, normal;
473  double angle, x, y;
474 
475  cross_product(v1, v2, &normal);
476  normalize(&normal);
477  cross_product(&normal, v1, &v3);
478 
479  x = dot_product(v1, v2);
480  y = dot_product(v2, &v3);
481 
482  angle = atan2(y, x);
483  return angle;
484 }
485 
489 static void normalize2d(POINT2D *p)
490 {
491  double d = sqrt(p->x*p->x + p->y*p->y);
492  if (FP_IS_ZERO(d))
493  {
494  p->x = p->y = 0.0;
495  return;
496  }
497  p->x = p->x / d;
498  p->y = p->y / d;
499  return;
500 }
501 
506 void unit_normal(const POINT3D *P1, const POINT3D *P2, POINT3D *normal)
507 {
508  double p_dot = dot_product(P1, P2);
509  POINT3D P3;
510 
511  /* If edge is really large, calculate a narrower equivalent angle A1/A3. */
512  if ( p_dot < 0 )
513  {
514  vector_sum(P1, P2, &P3);
515  normalize(&P3);
516  }
517  /* If edge is narrow, calculate a wider equivalent angle A1/A3. */
518  else if ( p_dot > 0.95 )
519  {
520  vector_difference(P2, P1, &P3);
521  normalize(&P3);
522  }
523  /* Just keep the current angle in A1/A3. */
524  else
525  {
526  P3 = *P2;
527  }
528 
529  /* Normals to the A-plane and B-plane */
530  cross_product(P1, &P3, normal);
531  normalize(normal);
532 }
533 
538 void vector_rotate(const POINT3D* v1, const POINT3D* v2, double angle, POINT3D* n)
539 {
540  POINT3D u;
541  double cos_a = cos(angle);
542  double sin_a = sin(angle);
543  double uxuy, uyuz, uxuz;
544  double ux2, uy2, uz2;
545  double rxx, rxy, rxz, ryx, ryy, ryz, rzx, rzy, rzz;
546 
547  /* Need a unit vector normal to rotate around */
548  unit_normal(v1, v2, &u);
549 
550  uxuy = u.x * u.y;
551  uxuz = u.x * u.z;
552  uyuz = u.y * u.z;
553 
554  ux2 = u.x * u.x;
555  uy2 = u.y * u.y;
556  uz2 = u.z * u.z;
557 
558  rxx = cos_a + ux2 * (1 - cos_a);
559  rxy = uxuy * (1 - cos_a) - u.z * sin_a;
560  rxz = uxuz * (1 - cos_a) + u.y * sin_a;
561 
562  ryx = uxuy * (1 - cos_a) + u.z * sin_a;
563  ryy = cos_a + uy2 * (1 - cos_a);
564  ryz = uyuz * (1 - cos_a) - u.x * sin_a;
565 
566  rzx = uxuz * (1 - cos_a) - u.y * sin_a;
567  rzy = uyuz * (1 - cos_a) + u.x * sin_a;
568  rzz = cos_a + uz2 * (1 - cos_a);
569 
570  n->x = rxx * v1->x + rxy * v1->y + rxz * v1->z;
571  n->y = ryx * v1->x + ryy * v1->y + ryz * v1->z;
572  n->z = rzx * v1->x + rzy * v1->y + rzz * v1->z;
573 
574  normalize(n);
575 }
576 
581 {
582  double d = sqrt(p->x*p->x + p->y*p->y + p->z*p->z);
583  if (FP_IS_ZERO(d))
584  {
585  p->x = p->y = p->z = 0.0;
586  return;
587  }
588  p->x = p->x / d;
589  p->y = p->y / d;
590  p->z = p->z / d;
591  return;
592 }
593 
594 
600 {
601  double lon_qpp = (q->lon + p->lon) / -2.0;
602  double lon_qmp = (q->lon - p->lon) / 2.0;
603  double sin_p_lat_minus_q_lat = sin(p->lat-q->lat);
604  double sin_p_lat_plus_q_lat = sin(p->lat+q->lat);
605  double sin_lon_qpp = sin(lon_qpp);
606  double sin_lon_qmp = sin(lon_qmp);
607  double cos_lon_qpp = cos(lon_qpp);
608  double cos_lon_qmp = cos(lon_qmp);
609  a->x = sin_p_lat_minus_q_lat * sin_lon_qpp * cos_lon_qmp -
610  sin_p_lat_plus_q_lat * cos_lon_qpp * sin_lon_qmp;
611  a->y = sin_p_lat_minus_q_lat * cos_lon_qpp * cos_lon_qmp +
612  sin_p_lat_plus_q_lat * sin_lon_qpp * sin_lon_qmp;
613  a->z = cos(p->lat) * cos(q->lat) * sin(q->lon-p->lon);
614 }
615 
616 void x_to_z(POINT3D *p)
617 {
618  double tmp = p->z;
619  p->z = p->x;
620  p->x = tmp;
621 }
622 
623 void y_to_z(POINT3D *p)
624 {
625  double tmp = p->z;
626  p->z = p->y;
627  p->y = tmp;
628 }
629 
630 
632 {
633  double sign_s = SIGNUM(s->lon);
634  double sign_e = SIGNUM(e->lon);
635  double ss = fabs(s->lon);
636  double ee = fabs(e->lon);
637  if ( sign_s == sign_e )
638  {
639  return LW_FALSE;
640  }
641  else
642  {
643  double dl = ss + ee;
644  if ( dl < M_PI )
645  return LW_FALSE;
646  else if ( FP_EQUALS(dl, M_PI) )
647  return LW_FALSE;
648  else
649  return LW_TRUE;
650  }
651 }
652 
658 static int
660 {
661  POINT3D normal, pt;
662  double w;
663  /* Normal to the plane defined by e */
664  robust_cross_product(&(e->start), &(e->end), &normal);
665  normalize(&normal);
666  geog2cart(p, &pt);
667  /* We expect the dot product of with normal with any vector in the plane to be zero */
668  w = dot_product(&normal, &pt);
669  LWDEBUGF(4,"dot product %.9g",w);
670  if ( FP_IS_ZERO(w) )
671  {
672  LWDEBUG(4, "point is on plane (dot product is zero)");
673  return 0;
674  }
675 
676  if ( w < 0 )
677  return -1;
678  else
679  return 1;
680 }
681 
685 static double
687 {
688  POINT3D normal1, normal2;
689  robust_cross_product(b, a, &normal1);
690  robust_cross_product(b, c, &normal2);
691  normalize(&normal1);
692  normalize(&normal2);
693  return sphere_distance_cartesian(&normal1, &normal2);
694 }
695 
705 static double
707 {
708  double angle_a, angle_b, angle_c;
709  double area_radians = 0.0;
710  int side;
711  GEOGRAPHIC_EDGE e;
712 
713  angle_a = sphere_angle(b,a,c);
714  angle_b = sphere_angle(a,b,c);
715  angle_c = sphere_angle(b,c,a);
716 
717  area_radians = angle_a + angle_b + angle_c - M_PI;
718 
719  /* What's the direction of the B/C edge? */
720  e.start = *a;
721  e.end = *b;
722  side = edge_point_side(&e, c);
723 
724  /* Co-linear points implies no area */
725  if ( side == 0 )
726  return 0.0;
727 
728  /* Add the sign to the area */
729  return side * area_radians;
730 }
731 
732 
733 
741 {
742  int side = edge_point_side(e, p);
743  if ( side == 0 )
744  return LW_TRUE;
745 
746  return LW_FALSE;
747 }
748 
754 {
755  POINT3D vcp, vs, ve, vp;
756  double vs_dot_vcp, vp_dot_vcp;
757  geog2cart(&(e->start), &vs);
758  geog2cart(&(e->end), &ve);
759  /* Antipodal case, everything is inside. */
760  if ( vs.x == -1.0 * ve.x && vs.y == -1.0 * ve.y && vs.z == -1.0 * ve.z )
761  return LW_TRUE;
762  geog2cart(p, &vp);
763  /* The normalized sum bisects the angle between start and end. */
764  vector_sum(&vs, &ve, &vcp);
765  normalize(&vcp);
766  /* The projection of start onto the center defines the minimum similarity */
767  vs_dot_vcp = dot_product(&vs, &vcp);
768  LWDEBUGF(4,"vs_dot_vcp %.19g",vs_dot_vcp);
769  /* The projection of candidate p onto the center */
770  vp_dot_vcp = dot_product(&vp, &vcp);
771  LWDEBUGF(4,"vp_dot_vcp %.19g",vp_dot_vcp);
772  /* If p is more similar than start then p is inside the cone */
773  LWDEBUGF(4,"fabs(vp_dot_vcp - vs_dot_vcp) %.39g",fabs(vp_dot_vcp - vs_dot_vcp));
774 
775  /*
776  ** We want to test that vp_dot_vcp is >= vs_dot_vcp but there are
777  ** numerical stability issues for values that are very very nearly
778  ** equal. Unfortunately there are also values of vp_dot_vcp that are legitimately
779  ** very close to but still less than vs_dot_vcp which we also need to catch.
780  ** The tolerance of 10-17 seems to do the trick on 32-bit and 64-bit architectures,
781  ** for the test cases here.
782  ** However, tuning the tolerance value feels like a dangerous hack.
783  ** Fundamentally, the problem is that this test is so sensitive.
784  */
785 
786  /* 1.1102230246251565404236316680908203125e-16 */
787 
788  if ( vp_dot_vcp > vs_dot_vcp || fabs(vp_dot_vcp - vs_dot_vcp) < 2e-16 )
789  {
790  LWDEBUG(4, "point is in cone");
791  return LW_TRUE;
792  }
793  LWDEBUG(4, "point is not in cone");
794  return LW_FALSE;
795 }
796 
801 {
802  GEOGRAPHIC_EDGE g;
804  double slon = fabs((e->start).lon) + fabs((e->end).lon);
805  double dlon = fabs(fabs((e->start).lon) - fabs((e->end).lon));
806  double slat = (e->start).lat + (e->end).lat;
807 
808  LWDEBUGF(4, "e.start == GPOINT(%.6g %.6g) ", (e->start).lat, (e->start).lon);
809  LWDEBUGF(4, "e.end == GPOINT(%.6g %.6g) ", (e->end).lat, (e->end).lon);
810  LWDEBUGF(4, "p == GPOINT(%.6g %.6g) ", p->lat, p->lon);
811 
812  /* Copy values into working registers */
813  g = *e;
814  q = *p;
815 
816  /* Vertical plane, we need to do this calculation in latitude */
817  if ( FP_EQUALS( g.start.lon, g.end.lon ) )
818  {
819  LWDEBUG(4, "vertical plane, we need to do this calculation in latitude");
820  /* Supposed to be co-planar... */
821  if ( ! FP_EQUALS( q.lon, g.start.lon ) )
822  return LW_FALSE;
823 
824  if ( ( g.start.lat <= q.lat && q.lat <= g.end.lat ) ||
825  ( g.end.lat <= q.lat && q.lat <= g.start.lat ) )
826  {
827  return LW_TRUE;
828  }
829  else
830  {
831  return LW_FALSE;
832  }
833  }
834 
835  /* Over the pole, we need normalize latitude and do this calculation in latitude */
836  if ( FP_EQUALS( slon, M_PI ) && ( SIGNUM(g.start.lon) != SIGNUM(g.end.lon) || FP_EQUALS(dlon, M_PI) ) )
837  {
838  LWDEBUG(4, "over the pole...");
839  /* Antipodal, everything (or nothing?) is inside */
840  if ( FP_EQUALS( slat, 0.0 ) )
841  return LW_TRUE;
842 
843  /* Point *is* the north pole */
844  if ( slat > 0.0 && FP_EQUALS(q.lat, M_PI_2 ) )
845  return LW_TRUE;
846 
847  /* Point *is* the south pole */
848  if ( slat < 0.0 && FP_EQUALS(q.lat, -1.0 * M_PI_2) )
849  return LW_TRUE;
850 
851  LWDEBUG(4, "coplanar?...");
852 
853  /* Supposed to be co-planar... */
854  if ( ! FP_EQUALS( q.lon, g.start.lon ) )
855  return LW_FALSE;
856 
857  LWDEBUG(4, "north or south?...");
858 
859  /* Over north pole, test based on south pole */
860  if ( slat > 0.0 )
861  {
862  LWDEBUG(4, "over the north pole...");
863  if ( q.lat > FP_MIN(g.start.lat, g.end.lat) )
864  return LW_TRUE;
865  else
866  return LW_FALSE;
867  }
868  else
869  /* Over south pole, test based on north pole */
870  {
871  LWDEBUG(4, "over the south pole...");
872  if ( q.lat < FP_MAX(g.start.lat, g.end.lat) )
873  return LW_TRUE;
874  else
875  return LW_FALSE;
876  }
877  }
878 
879  /* Dateline crossing, flip everything to the opposite hemisphere */
880  else if ( slon > M_PI && ( SIGNUM(g.start.lon) != SIGNUM(g.end.lon) ) )
881  {
882  LWDEBUG(4, "crosses dateline, flip longitudes...");
883  if ( g.start.lon > 0.0 )
884  g.start.lon -= M_PI;
885  else
886  g.start.lon += M_PI;
887  if ( g.end.lon > 0.0 )
888  g.end.lon -= M_PI;
889  else
890  g.end.lon += M_PI;
891 
892  if ( q.lon > 0.0 )
893  q.lon -= M_PI;
894  else
895  q.lon += M_PI;
896  }
897 
898  if ( ( g.start.lon <= q.lon && q.lon <= g.end.lon ) ||
899  ( g.end.lon <= q.lon && q.lon <= g.start.lon ) )
900  {
901  LWDEBUG(4, "true, this edge contains point");
902  return LW_TRUE;
903  }
904 
905  LWDEBUG(4, "false, this edge does not contain point");
906  return LW_FALSE;
907 }
908 
909 
914 {
915  double d_lon = e->lon - s->lon;
916  double cos_d_lon = cos(d_lon);
917  double cos_lat_e = cos(e->lat);
918  double sin_lat_e = sin(e->lat);
919  double cos_lat_s = cos(s->lat);
920  double sin_lat_s = sin(s->lat);
921 
922  double a1 = POW2(cos_lat_e * sin(d_lon));
923  double a2 = POW2(cos_lat_s * sin_lat_e - sin_lat_s * cos_lat_e * cos_d_lon);
924  double a = sqrt(a1 + a2);
925  double b = sin_lat_s * sin_lat_e + cos_lat_s * cos_lat_e * cos_d_lon;
926  return atan2(a, b);
927 }
928 
932 double sphere_distance_cartesian(const POINT3D *s, const POINT3D *e)
933 {
934  return acos(FP_MIN(1.0, dot_product(s, e)));
935 }
936 
940 double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, double d)
941 {
942  double heading = 0.0;
943  double f;
944 
945  /* Starting from the poles? Special case. */
946  if ( FP_IS_ZERO(cos(s->lat)) )
947  return (s->lat > 0.0) ? M_PI : 0.0;
948 
949  f = (sin(e->lat) - sin(s->lat) * cos(d)) / (sin(d) * cos(s->lat));
950  if ( FP_EQUALS(f, 1.0) )
951  heading = 0.0;
952  else if ( FP_EQUALS(f, -1.0) )
953  heading = M_PI;
954  else if ( fabs(f) > 1.0 )
955  {
956  LWDEBUGF(4, "f = %g", f);
957  heading = acos(f);
958  }
959  else
960  heading = acos(f);
961 
962  if ( sin(e->lon - s->lon) < 0.0 )
963  heading = -1 * heading;
964 
965  return heading;
966 }
967 
968 #if 0 /* unused */
969 
980 static double sphere_excess(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c)
981 {
982  double a_dist = sphere_distance(b, c);
983  double b_dist = sphere_distance(c, a);
984  double c_dist = sphere_distance(a, b);
985  double hca = sphere_direction(c, a, b_dist);
986  double hcb = sphere_direction(c, b, a_dist);
987  double sign = SIGNUM(hcb-hca);
988  double ss = (a_dist + b_dist + c_dist) / 2.0;
989  double E = tan(ss/2.0)*tan((ss-a_dist)/2.0)*tan((ss-b_dist)/2.0)*tan((ss-c_dist)/2.0);
990  return 4.0 * atan(sqrt(fabs(E))) * sign;
991 }
992 #endif
993 
994 
1000 {
1001  if ( edge_point_in_cone(e, p) && edge_point_on_plane(e, p) )
1002  /* if ( edge_contains_coplanar_point(e, p) && edge_point_on_plane(e, p) ) */
1003  {
1004  LWDEBUG(4, "point is on edge");
1005  return LW_TRUE;
1006  }
1007  LWDEBUG(4, "point is not on edge");
1008  return LW_FALSE;
1009 }
1010 
1014 double z_to_latitude(double z, int top)
1015 {
1016  double sign = SIGNUM(z);
1017  double tlat = acos(z);
1018  LWDEBUGF(4, "inputs: z(%.8g) sign(%.8g) tlat(%.8g)", z, sign, tlat);
1019  if (FP_IS_ZERO(z))
1020  {
1021  if (top) return M_PI_2;
1022  else return -1.0 * M_PI_2;
1023  }
1024  if (fabs(tlat) > M_PI_2 )
1025  {
1026  tlat = sign * (M_PI - fabs(tlat));
1027  }
1028  else
1029  {
1030  tlat = sign * tlat;
1031  }
1032  LWDEBUGF(4, "output: tlat(%.8g)", tlat);
1033  return tlat;
1034 }
1035 
1041 int clairaut_cartesian(const POINT3D *start, const POINT3D *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
1042 {
1043  POINT3D t1, t2;
1044  GEOGRAPHIC_POINT vN1, vN2;
1045  LWDEBUG(4,"entering function");
1046  unit_normal(start, end, &t1);
1047  unit_normal(end, start, &t2);
1048  LWDEBUGF(4, "unit normal t1 == POINT(%.8g %.8g %.8g)", t1.x, t1.y, t1.z);
1049  LWDEBUGF(4, "unit normal t2 == POINT(%.8g %.8g %.8g)", t2.x, t2.y, t2.z);
1050  cart2geog(&t1, &vN1);
1051  cart2geog(&t2, &vN2);
1052  g_top->lat = z_to_latitude(t1.z,LW_TRUE);
1053  g_top->lon = vN2.lon;
1054  g_bottom->lat = z_to_latitude(t2.z,LW_FALSE);
1055  g_bottom->lon = vN1.lon;
1056  LWDEBUGF(4, "clairaut top == GPOINT(%.6g %.6g)", g_top->lat, g_top->lon);
1057  LWDEBUGF(4, "clairaut bottom == GPOINT(%.6g %.6g)", g_bottom->lat, g_bottom->lon);
1058  return LW_SUCCESS;
1059 }
1060 
1067 {
1068  POINT3D t1, t2;
1069  GEOGRAPHIC_POINT vN1, vN2;
1070  LWDEBUG(4,"entering function");
1071  robust_cross_product(start, end, &t1);
1072  normalize(&t1);
1073  robust_cross_product(end, start, &t2);
1074  normalize(&t2);
1075  LWDEBUGF(4, "unit normal t1 == POINT(%.8g %.8g %.8g)", t1.x, t1.y, t1.z);
1076  LWDEBUGF(4, "unit normal t2 == POINT(%.8g %.8g %.8g)", t2.x, t2.y, t2.z);
1077  cart2geog(&t1, &vN1);
1078  cart2geog(&t2, &vN2);
1079  g_top->lat = z_to_latitude(t1.z,LW_TRUE);
1080  g_top->lon = vN2.lon;
1081  g_bottom->lat = z_to_latitude(t2.z,LW_FALSE);
1082  g_bottom->lon = vN1.lon;
1083  LWDEBUGF(4, "clairaut top == GPOINT(%.6g %.6g)", g_top->lat, g_top->lon);
1084  LWDEBUGF(4, "clairaut bottom == GPOINT(%.6g %.6g)", g_bottom->lat, g_bottom->lon);
1085  return LW_SUCCESS;
1086 }
1087 
1093 {
1094  POINT3D ea, eb, v;
1095  LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", e1->start.lat, e1->start.lon, e1->end.lat, e1->end.lon);
1096  LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", e2->start.lat, e2->start.lon, e2->end.lat, e2->end.lon);
1097 
1098  LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e1->start.lon), rad2deg(e1->start.lat), rad2deg(e1->end.lon), rad2deg(e1->end.lat));
1099  LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e2->start.lon), rad2deg(e2->start.lat), rad2deg(e2->end.lon), rad2deg(e2->end.lat));
1100 
1101  if ( geographic_point_equals(&(e1->start), &(e2->start)) )
1102  {
1103  *g = e1->start;
1104  return LW_TRUE;
1105  }
1106  if ( geographic_point_equals(&(e1->end), &(e2->end)) )
1107  {
1108  *g = e1->end;
1109  return LW_TRUE;
1110  }
1111  if ( geographic_point_equals(&(e1->end), &(e2->start)) )
1112  {
1113  *g = e1->end;
1114  return LW_TRUE;
1115  }
1116  if ( geographic_point_equals(&(e1->start), &(e2->end)) )
1117  {
1118  *g = e1->start;
1119  return LW_TRUE;
1120  }
1121 
1122  robust_cross_product(&(e1->start), &(e1->end), &ea);
1123  normalize(&ea);
1124  robust_cross_product(&(e2->start), &(e2->end), &eb);
1125  normalize(&eb);
1126  LWDEBUGF(4, "e1 cross product == POINT(%.12g %.12g %.12g)", ea.x, ea.y, ea.z);
1127  LWDEBUGF(4, "e2 cross product == POINT(%.12g %.12g %.12g)", eb.x, eb.y, eb.z);
1128  LWDEBUGF(4, "fabs(dot_product(ea, eb)) == %.14g", fabs(dot_product(&ea, &eb)));
1129  if ( FP_EQUALS(fabs(dot_product(&ea, &eb)), 1.0) )
1130  {
1131  LWDEBUGF(4, "parallel edges found! dot_product = %.12g", dot_product(&ea, &eb));
1132  /* Parallel (maybe equal) edges! */
1133  /* Hack alert, only returning ONE end of the edge right now, most do better later. */
1134  /* Hack alart #2, returning a value of 2 to indicate a co-linear crossing event. */
1135  if ( edge_contains_point(e1, &(e2->start)) )
1136  {
1137  *g = e2->start;
1138  return 2;
1139  }
1140  if ( edge_contains_point(e1, &(e2->end)) )
1141  {
1142  *g = e2->end;
1143  return 2;
1144  }
1145  if ( edge_contains_point(e2, &(e1->start)) )
1146  {
1147  *g = e1->start;
1148  return 2;
1149  }
1150  if ( edge_contains_point(e2, &(e1->end)) )
1151  {
1152  *g = e1->end;
1153  return 2;
1154  }
1155  }
1156  unit_normal(&ea, &eb, &v);
1157  LWDEBUGF(4, "v == POINT(%.12g %.12g %.12g)", v.x, v.y, v.z);
1158  g->lat = atan2(v.z, sqrt(v.x * v.x + v.y * v.y));
1159  g->lon = atan2(v.y, v.x);
1160  LWDEBUGF(4, "g == GPOINT(%.12g %.12g)", g->lat, g->lon);
1161  LWDEBUGF(4, "g == POINT(%.12g %.12g)", rad2deg(g->lon), rad2deg(g->lat));
1162  if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
1163  {
1164  return LW_TRUE;
1165  }
1166  else
1167  {
1168  LWDEBUG(4, "flipping point to other side of sphere");
1169  g->lat = -1.0 * g->lat;
1170  g->lon = g->lon + M_PI;
1171  if ( g->lon > M_PI )
1172  {
1173  g->lon = -1.0 * (2.0 * M_PI - g->lon);
1174  }
1175  if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
1176  {
1177  return LW_TRUE;
1178  }
1179  }
1180  return LW_FALSE;
1181 }
1182 
1184 {
1185  double d1 = 1000000000.0, d2, d3, d_nearest;
1186  POINT3D n, p, k;
1187  GEOGRAPHIC_POINT gk, g_nearest;
1188 
1189  /* Zero length edge, */
1190  if ( geographic_point_equals(&(e->start), &(e->end)) )
1191  {
1192  *closest = e->start;
1193  return sphere_distance(&(e->start), gp);
1194  }
1195 
1196  robust_cross_product(&(e->start), &(e->end), &n);
1197  normalize(&n);
1198  geog2cart(gp, &p);
1199  vector_scale(&n, dot_product(&p, &n));
1200  vector_difference(&p, &n, &k);
1201  normalize(&k);
1202  cart2geog(&k, &gk);
1203  if ( edge_contains_point(e, &gk) )
1204  {
1205  d1 = sphere_distance(gp, &gk);
1206  }
1207  d2 = sphere_distance(gp, &(e->start));
1208  d3 = sphere_distance(gp, &(e->end));
1209 
1210  d_nearest = d1;
1211  g_nearest = gk;
1212 
1213  if ( d2 < d_nearest )
1214  {
1215  d_nearest = d2;
1216  g_nearest = e->start;
1217  }
1218  if ( d3 < d_nearest )
1219  {
1220  d_nearest = d3;
1221  g_nearest = e->end;
1222  }
1223  if (closest)
1224  *closest = g_nearest;
1225 
1226  return d_nearest;
1227 }
1228 
1235 {
1236  double d;
1237  GEOGRAPHIC_POINT gcp1s, gcp1e, gcp2s, gcp2e, c1, c2;
1238  double d1s = edge_distance_to_point(e1, &(e2->start), &gcp1s);
1239  double d1e = edge_distance_to_point(e1, &(e2->end), &gcp1e);
1240  double d2s = edge_distance_to_point(e2, &(e1->start), &gcp2s);
1241  double d2e = edge_distance_to_point(e2, &(e1->end), &gcp2e);
1242 
1243  d = d1s;
1244  c1 = gcp1s;
1245  c2 = e2->start;
1246 
1247  if ( d1e < d )
1248  {
1249  d = d1e;
1250  c1 = gcp1e;
1251  c2 = e2->end;
1252  }
1253 
1254  if ( d2s < d )
1255  {
1256  d = d2s;
1257  c1 = e1->start;
1258  c2 = gcp2s;
1259  }
1260 
1261  if ( d2e < d )
1262  {
1263  d = d2e;
1264  c1 = e1->end;
1265  c2 = gcp2e;
1266  }
1267 
1268  if ( closest1 ) *closest1 = c1;
1269  if ( closest2 ) *closest2 = c2;
1270 
1271  return d;
1272 }
1273 
1274 
1279 int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, GEOGRAPHIC_POINT *n)
1280 {
1281  double d = distance;
1282  double lat1 = r->lat;
1283  double lon1 = r->lon;
1284  double lat2, lon2;
1285 
1286  lat2 = asin(sin(lat1)*cos(d) + cos(lat1)*sin(d)*cos(azimuth));
1287 
1288  /* If we're going straight up or straight down, we don't need to calculate the longitude */
1289  /* TODO: this isn't quite true, what if we're going over the pole? */
1290  if ( FP_EQUALS(azimuth, M_PI) || FP_EQUALS(azimuth, 0.0) )
1291  {
1292  lon2 = r->lon;
1293  }
1294  else
1295  {
1296  lon2 = lon1 + atan2(sin(azimuth)*sin(d)*cos(lat1), cos(d)-sin(lat1)*sin(lat2));
1297  }
1298 
1299  if ( isnan(lat2) || isnan(lon2) )
1300  return LW_FAILURE;
1301 
1302  n->lat = lat2;
1303  n->lon = lon2;
1304 
1305  return LW_SUCCESS;
1306 }
1307 
1308 
1310 {
1311  int steps = 1000000;
1312  int i;
1313  double dx, dy, dz;
1314  double distance = sphere_distance(&(e->start), &(e->end));
1315  POINT3D pn, p, start, end;
1316 
1317  /* Edge is zero length, just return the naive box */
1318  if ( FP_IS_ZERO(distance) )
1319  {
1320  LWDEBUG(4, "edge is zero length. returning");
1321  geog2cart(&(e->start), &start);
1322  geog2cart(&(e->end), &end);
1323  gbox_init_point3d(&start, gbox);
1324  gbox_merge_point3d(&end, gbox);
1325  return LW_SUCCESS;
1326  }
1327 
1328  /* Edge is antipodal (one point on each side of the globe),
1329  set the box to contain the whole world and return */
1330  if ( FP_EQUALS(distance, M_PI) )
1331  {
1332  LWDEBUG(4, "edge is antipodal. setting to maximum size box, and returning");
1333  gbox->xmin = gbox->ymin = gbox->zmin = -1.0;
1334  gbox->xmax = gbox->ymax = gbox->zmax = 1.0;
1335  return LW_SUCCESS;
1336  }
1337 
1338  /* Walk along the chord between start and end incrementally,
1339  normalizing at each step. */
1340  geog2cart(&(e->start), &start);
1341  geog2cart(&(e->end), &end);
1342  dx = (end.x - start.x)/steps;
1343  dy = (end.y - start.y)/steps;
1344  dz = (end.z - start.z)/steps;
1345  p = start;
1346  gbox->xmin = gbox->xmax = p.x;
1347  gbox->ymin = gbox->ymax = p.y;
1348  gbox->zmin = gbox->zmax = p.z;
1349  for ( i = 0; i < steps; i++ )
1350  {
1351  p.x += dx;
1352  p.y += dy;
1353  p.z += dz;
1354  pn = p;
1355  normalize(&pn);
1356  gbox_merge_point3d(&pn, gbox);
1357  }
1358  return LW_SUCCESS;
1359 }
1360 
1373 int edge_calculate_gbox(const POINT3D *A1, const POINT3D *A2, GBOX *gbox)
1374 {
1375  POINT2D R1, R2, RX, O;
1376  POINT3D AN, A3;
1377  POINT3D X[6];
1378  int i, o_side;
1379 
1380  /* Initialize the box with the edge end points */
1381  gbox_init_point3d(A1, gbox);
1382  gbox_merge_point3d(A2, gbox);
1383 
1384  /* Zero length edge, just return! */
1385  if ( p3d_same(A1, A2) )
1386  return LW_SUCCESS;
1387 
1388  /* Error out on antipodal edge */
1389  if ( FP_EQUALS(A1->x, -1*A2->x) && FP_EQUALS(A1->y, -1*A2->y) && FP_EQUALS(A1->z, -1*A2->z) )
1390  {
1391  lwerror("Antipodal (180 degrees long) edge detected!");
1392  return LW_FAILURE;
1393  }
1394 
1395  /* Create A3, a vector in the plane of A1/A2, orthogonal to A1 */
1396  unit_normal(A1, A2, &AN);
1397  unit_normal(&AN, A1, &A3);
1398 
1399  /* Project A1 and A2 into the 2-space formed by the plane A1/A3 */
1400  R1.x = 1.0;
1401  R1.y = 0.0;
1402  R2.x = dot_product(A2, A1);
1403  R2.y = dot_product(A2, &A3);
1404 
1405  /* Initialize our 3-space axis points (x+, x-, y+, y-, z+, z-) */
1406  memset(X, 0, sizeof(POINT3D) * 6);
1407  X[0].x = X[2].y = X[4].z = 1.0;
1408  X[1].x = X[3].y = X[5].z = -1.0;
1409 
1410  /* Initialize a 2-space origin point. */
1411  O.x = O.y = 0.0;
1412  /* What side of the line joining R1/R2 is O? */
1413  o_side = lw_segment_side(&R1, &R2, &O);
1414 
1415  /* Add any extrema! */
1416  for ( i = 0; i < 6; i++ )
1417  {
1418  /* Convert 3-space axis points to 2-space unit vectors */
1419  RX.x = dot_product(&(X[i]), A1);
1420  RX.y = dot_product(&(X[i]), &A3);
1421  normalize2d(&RX);
1422 
1423  /* Any axis end on the side of R1/R2 opposite the origin */
1424  /* is an extreme point in the arc, so we add the 3-space */
1425  /* version of the point on R1/R2 to the gbox */
1426  if ( lw_segment_side(&R1, &R2, &RX) != o_side )
1427  {
1428  POINT3D Xn;
1429  Xn.x = RX.x * A1->x + RX.y * A3.x;
1430  Xn.y = RX.x * A1->y + RX.y * A3.y;
1431  Xn.z = RX.x * A1->z + RX.y * A3.z;
1432 
1433  gbox_merge_point3d(&Xn, gbox);
1434  }
1435  }
1436 
1437  return LW_SUCCESS;
1438 }
1439 
1440 void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside)
1441 {
1442  /* Make sure we have boxes */
1443  if ( poly->bbox )
1444  {
1445  gbox_pt_outside(poly->bbox, pt_outside);
1446  return;
1447  }
1448  else
1449  {
1450  GBOX gbox;
1451  lwgeom_calculate_gbox_geodetic((LWGEOM*)poly, &gbox);
1452  gbox_pt_outside(&gbox, pt_outside);
1453  return;
1454  }
1455 }
1456 
1461 void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
1462 {
1463  double grow = M_PI / 180.0 / 60.0; /* one arc-minute */
1464  int i;
1465  GBOX ge;
1466  POINT3D corners[8];
1467  POINT3D pt;
1468  GEOGRAPHIC_POINT g;
1469 
1470  while ( grow < M_PI )
1471  {
1472  /* Assign our box and expand it slightly. */
1473  ge = *gbox;
1474  if ( ge.xmin > -1 ) ge.xmin -= grow;
1475  if ( ge.ymin > -1 ) ge.ymin -= grow;
1476  if ( ge.zmin > -1 ) ge.zmin -= grow;
1477  if ( ge.xmax < 1 ) ge.xmax += grow;
1478  if ( ge.ymax < 1 ) ge.ymax += grow;
1479  if ( ge.zmax < 1 ) ge.zmax += grow;
1480 
1481  /* Build our eight corner points */
1482  corners[0].x = ge.xmin;
1483  corners[0].y = ge.ymin;
1484  corners[0].z = ge.zmin;
1485 
1486  corners[1].x = ge.xmin;
1487  corners[1].y = ge.ymax;
1488  corners[1].z = ge.zmin;
1489 
1490  corners[2].x = ge.xmin;
1491  corners[2].y = ge.ymin;
1492  corners[2].z = ge.zmax;
1493 
1494  corners[3].x = ge.xmax;
1495  corners[3].y = ge.ymin;
1496  corners[3].z = ge.zmin;
1497 
1498  corners[4].x = ge.xmax;
1499  corners[4].y = ge.ymax;
1500  corners[4].z = ge.zmin;
1501 
1502  corners[5].x = ge.xmax;
1503  corners[5].y = ge.ymin;
1504  corners[5].z = ge.zmax;
1505 
1506  corners[6].x = ge.xmin;
1507  corners[6].y = ge.ymax;
1508  corners[6].z = ge.zmax;
1509 
1510  corners[7].x = ge.xmax;
1511  corners[7].y = ge.ymax;
1512  corners[7].z = ge.zmax;
1513 
1514  LWDEBUG(4, "trying to use a box corner point...");
1515  for ( i = 0; i < 8; i++ )
1516  {
1517  normalize(&(corners[i]));
1518  LWDEBUGF(4, "testing corner %d: POINT(%.8g %.8g %.8g)", i, corners[i].x, corners[i].y, corners[i].z);
1519  if ( ! gbox_contains_point3d(gbox, &(corners[i])) )
1520  {
1521  LWDEBUGF(4, "corner %d is outside our gbox", i);
1522  pt = corners[i];
1523  normalize(&pt);
1524  cart2geog(&pt, &g);
1525  pt_outside->x = rad2deg(g.lon);
1526  pt_outside->y = rad2deg(g.lat);
1527  LWDEBUGF(4, "returning POINT(%.8g %.8g) as outside point", pt_outside->x, pt_outside->y);
1528  return;
1529  }
1530  }
1531 
1532  /* Try a wider growth to push the corners outside the original box. */
1533  grow *= 2.0;
1534  }
1535 
1536  /* This should never happen! */
1537  lwerror("BOOM! Could not generate outside point!");
1538  return;
1539 }
1540 
1541 
1543  const POINT3D *p1, const POINT3D *p2, /* 3-space points we are interpolating beween */
1544  const POINT4D *v1, const POINT4D *v2, /* real values and z/m values */
1545  double d, double max_seg_length, /* current segment length and segment limit */
1546  POINTARRAY *pa) /* write out results here */
1547 {
1548  /* Reached the terminal leaf in recursion. Add */
1549  /* the left-most point to the pointarray here */
1550  /* We recurse down the left side first, so outputs should */
1551  /* end up added to the array in order this way */
1552  if (d <= max_seg_length)
1553  {
1554  GEOGRAPHIC_POINT g;
1555  POINT4D p;
1556  cart2geog(p1, &g);
1557  p.x = rad2deg(g.lon);
1558  p.y = rad2deg(g.lat);
1559  p.z = v1->z;
1560  p.m = v1->m;
1561  return ptarray_append_point(pa, &p, LW_FALSE);
1562  }
1563  /* Find the mid-point and recurse on the left and then the right */
1564  else
1565  {
1566  /* Calculate mid-point */
1567  POINT3D mid;
1568  mid.x = (p1->x + p2->x) / 2.0;
1569  mid.y = (p1->y + p2->y) / 2.0;
1570  mid.z = (p1->z + p2->z) / 2.0;
1571  normalize(&mid);
1572 
1573  /* Calculate z/m mid-values */
1574  /* (ignore x/y, we get those from the 3-space calculations) */
1575  POINT4D midv;
1576  midv.z = (v1->z + v2->z) / 2.0;
1577  midv.m = (v1->m + v2->m) / 2.0;
1578  /* Recurse on the left first */
1579  ptarray_segmentize_sphere_edge_recursive(p1, &mid, v1, &midv, d/2.0, max_seg_length, pa);
1580  ptarray_segmentize_sphere_edge_recursive(&mid, p2, &midv, v2, d/2.0, max_seg_length, pa);
1581  return LW_SUCCESS;
1582  }
1583 }
1584 
1590 static POINTARRAY*
1591 ptarray_segmentize_sphere(const POINTARRAY *pa_in, double max_seg_length)
1592 {
1593  POINTARRAY *pa_out;
1594  int hasz = ptarray_has_z(pa_in);
1595  int hasm = ptarray_has_m(pa_in);
1596  POINT4D p1, p2;
1597  POINT3D q1, q2;
1598  GEOGRAPHIC_POINT g1, g2;
1599  int i;
1600 
1601  /* Just crap out on crazy input */
1602  if ( ! pa_in )
1603  lwerror("%s: null input pointarray", __func__);
1604  if ( max_seg_length <= 0.0 )
1605  lwerror("%s: maximum segment length must be positive", __func__);
1606 
1607  /* Empty starting array */
1608  pa_out = ptarray_construct_empty(hasz, hasm, pa_in->npoints);
1609 
1610  /* Simple loop per edge */
1611  for (i = 1; i < pa_in->npoints; i++)
1612  {
1613  getPoint4d_p(pa_in, i-1, &p1);
1614  getPoint4d_p(pa_in, i, &p2);
1615  geographic_point_init(p1.x, p1.y, &g1);
1616  geographic_point_init(p2.x, p2.y, &g2);
1617 
1618  /* Skip duplicate points (except in case of 2-point lines!) */
1619  if ((pa_in->npoints > 2) && p4d_same(&p1, &p2))
1620  continue;
1621 
1622  /* How long is this edge? */
1623  double d = sphere_distance(&g1, &g2);
1624 
1625  if (d > max_seg_length)
1626  {
1627  geog2cart(&g1, &q1);
1628  geog2cart(&g2, &q2);
1629  /* 3-d end points, XYZM end point, current edge size, min edge size */
1630  ptarray_segmentize_sphere_edge_recursive(&q1, &q2, &p1, &p2, d, max_seg_length, pa_out);
1631  }
1632  /* If we don't segmentize, we need to add first point manually */
1633  else
1634  {
1635  ptarray_append_point(pa_out, &p1, LW_TRUE);
1636  }
1637  }
1638  /* Always add the last point */
1639  ptarray_append_point(pa_out, &p2, LW_TRUE);
1640  return pa_out;
1641 }
1642 
1649 LWGEOM*
1650 lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length)
1651 {
1652  POINTARRAY *pa_out;
1653  LWLINE *lwline;
1654  LWPOLY *lwpoly_in, *lwpoly_out;
1655  LWCOLLECTION *lwcol_in, *lwcol_out;
1656  int i;
1657 
1658  /* Reflect NULL */
1659  if ( ! lwg_in )
1660  return NULL;
1661 
1662  /* Clone empty */
1663  if ( lwgeom_is_empty(lwg_in) )
1664  return lwgeom_clone(lwg_in);
1665 
1666  switch (lwg_in->type)
1667  {
1668  case MULTIPOINTTYPE:
1669  case POINTTYPE:
1670  return lwgeom_clone_deep(lwg_in);
1671  break;
1672  case LINETYPE:
1673  lwline = lwgeom_as_lwline(lwg_in);
1674  pa_out = ptarray_segmentize_sphere(lwline->points, max_seg_length);
1675  return lwline_as_lwgeom(lwline_construct(lwg_in->srid, NULL, pa_out));
1676  break;
1677  case POLYGONTYPE:
1678  lwpoly_in = lwgeom_as_lwpoly(lwg_in);
1679  lwpoly_out = lwpoly_construct_empty(lwg_in->srid, lwgeom_has_z(lwg_in), lwgeom_has_m(lwg_in));
1680  for ( i = 0; i < lwpoly_in->nrings; i++ )
1681  {
1682  pa_out = ptarray_segmentize_sphere(lwpoly_in->rings[i], max_seg_length);
1683  lwpoly_add_ring(lwpoly_out, pa_out);
1684  }
1685  return lwpoly_as_lwgeom(lwpoly_out);
1686  break;
1687  case MULTILINETYPE:
1688  case MULTIPOLYGONTYPE:
1689  case COLLECTIONTYPE:
1690  lwcol_in = lwgeom_as_lwcollection(lwg_in);
1691  lwcol_out = lwcollection_construct_empty(lwg_in->type, lwg_in->srid, lwgeom_has_z(lwg_in), lwgeom_has_m(lwg_in));
1692  for ( i = 0; i < lwcol_in->ngeoms; i++ )
1693  {
1694  lwcollection_add_lwgeom(lwcol_out, lwgeom_segmentize_sphere(lwcol_in->geoms[i], max_seg_length));
1695  }
1696  return lwcollection_as_lwgeom(lwcol_out);
1697  break;
1698  default:
1699  lwerror("lwgeom_segmentize_sphere: unsupported input geometry type: %d - %s",
1700  lwg_in->type, lwtype_name(lwg_in->type));
1701  break;
1702  }
1703 
1704  lwerror("lwgeom_segmentize_sphere got to the end of the function, should not happen");
1705  return NULL;
1706 }
1707 
1708 
1713 double
1715 {
1716  int i;
1717  const POINT2D *p;
1718  GEOGRAPHIC_POINT a, b, c;
1719  double area = 0.0;
1720 
1721  /* Return zero on nonsensical inputs */
1722  if ( ! pa || pa->npoints < 4 )
1723  return 0.0;
1724 
1725  p = getPoint2d_cp(pa, 0);
1726  geographic_point_init(p->x, p->y, &a);
1727  p = getPoint2d_cp(pa, 1);
1728  geographic_point_init(p->x, p->y, &b);
1729 
1730  for ( i = 2; i < pa->npoints-1; i++ )
1731  {
1732  p = getPoint2d_cp(pa, i);
1733  geographic_point_init(p->x, p->y, &c);
1734  area += sphere_signed_area(&a, &b, &c);
1735  b = c;
1736  }
1737 
1738  return fabs(area);
1739 }
1740 
1741 
1742 static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY *pa2, const SPHEROID *s, double tolerance, int check_intersection)
1743 {
1744  GEOGRAPHIC_EDGE e1, e2;
1745  GEOGRAPHIC_POINT g1, g2;
1746  GEOGRAPHIC_POINT nearest1, nearest2;
1747  POINT3D A1, A2, B1, B2;
1748  const POINT2D *p;
1749  double distance;
1750  int i, j;
1751  int use_sphere = (s->a == s->b ? 1 : 0);
1752 
1753  /* Make result really big, so that everything will be smaller than it */
1754  distance = FLT_MAX;
1755 
1756  /* Empty point arrays? Return negative */
1757  if ( pa1->npoints == 0 || pa2->npoints == 0 )
1758  return -1.0;
1759 
1760  /* Handle point/point case here */
1761  if ( pa1->npoints == 1 && pa2->npoints == 1 )
1762  {
1763  p = getPoint2d_cp(pa1, 0);
1764  geographic_point_init(p->x, p->y, &g1);
1765  p = getPoint2d_cp(pa2, 0);
1766  geographic_point_init(p->x, p->y, &g2);
1767  /* Sphere special case, axes equal */
1768  distance = s->radius * sphere_distance(&g1, &g2);
1769  if ( use_sphere )
1770  return distance;
1771  /* Below tolerance, actual distance isn't of interest */
1772  else if ( distance < 0.95 * tolerance )
1773  return distance;
1774  /* Close or greater than tolerance, get the real answer to be sure */
1775  else
1776  return spheroid_distance(&g1, &g2, s);
1777  }
1778 
1779  /* Handle point/line case here */
1780  if ( pa1->npoints == 1 || pa2->npoints == 1 )
1781  {
1782  /* Handle one/many case here */
1783  int i;
1784  const POINTARRAY *pa_one;
1785  const POINTARRAY *pa_many;
1786 
1787  if ( pa1->npoints == 1 )
1788  {
1789  pa_one = pa1;
1790  pa_many = pa2;
1791  }
1792  else
1793  {
1794  pa_one = pa2;
1795  pa_many = pa1;
1796  }
1797 
1798  /* Initialize our point */
1799  p = getPoint2d_cp(pa_one, 0);
1800  geographic_point_init(p->x, p->y, &g1);
1801 
1802  /* Initialize start of line */
1803  p = getPoint2d_cp(pa_many, 0);
1804  geographic_point_init(p->x, p->y, &(e1.start));
1805 
1806  /* Iterate through the edges in our line */
1807  for ( i = 1; i < pa_many->npoints; i++ )
1808  {
1809  double d;
1810  p = getPoint2d_cp(pa_many, i);
1811  geographic_point_init(p->x, p->y, &(e1.end));
1812  /* Get the spherical distance between point and edge */
1813  d = s->radius * edge_distance_to_point(&e1, &g1, &g2);
1814  /* New shortest distance! Record this distance / location */
1815  if ( d < distance )
1816  {
1817  distance = d;
1818  nearest2 = g2;
1819  }
1820  /* We've gotten closer than the tolerance... */
1821  if ( d < tolerance )
1822  {
1823  /* Working on a sphere? The answer is correct, return */
1824  if ( use_sphere )
1825  {
1826  return d;
1827  }
1828  /* Far enough past the tolerance that the spheroid calculation won't change things */
1829  else if ( d < tolerance * 0.95 )
1830  {
1831  return d;
1832  }
1833  /* On a spheroid and near the tolerance? Confirm that we are *actually* closer than tolerance */
1834  else
1835  {
1836  d = spheroid_distance(&g1, &nearest2, s);
1837  /* Yes, closer than tolerance, return! */
1838  if ( d < tolerance )
1839  return d;
1840  }
1841  }
1842  e1.start = e1.end;
1843  }
1844 
1845  /* On sphere, return answer */
1846  if ( use_sphere )
1847  return distance;
1848  /* On spheroid, calculate final answer based on closest approach */
1849  else
1850  return spheroid_distance(&g1, &nearest2, s);
1851 
1852  }
1853 
1854  /* Initialize start of line 1 */
1855  p = getPoint2d_cp(pa1, 0);
1856  geographic_point_init(p->x, p->y, &(e1.start));
1857  geog2cart(&(e1.start), &A1);
1858 
1859 
1860  /* Handle line/line case */
1861  for ( i = 1; i < pa1->npoints; i++ )
1862  {
1863  p = getPoint2d_cp(pa1, i);
1864  geographic_point_init(p->x, p->y, &(e1.end));
1865  geog2cart(&(e1.end), &A2);
1866 
1867  /* Initialize start of line 2 */
1868  p = getPoint2d_cp(pa2, 0);
1869  geographic_point_init(p->x, p->y, &(e2.start));
1870  geog2cart(&(e2.start), &B1);
1871 
1872  for ( j = 1; j < pa2->npoints; j++ )
1873  {
1874  double d;
1875 
1876  p = getPoint2d_cp(pa2, j);
1877  geographic_point_init(p->x, p->y, &(e2.end));
1878  geog2cart(&(e2.end), &B2);
1879 
1880  LWDEBUGF(4, "e1.start == GPOINT(%.6g %.6g) ", e1.start.lat, e1.start.lon);
1881  LWDEBUGF(4, "e1.end == GPOINT(%.6g %.6g) ", e1.end.lat, e1.end.lon);
1882  LWDEBUGF(4, "e2.start == GPOINT(%.6g %.6g) ", e2.start.lat, e2.start.lon);
1883  LWDEBUGF(4, "e2.end == GPOINT(%.6g %.6g) ", e2.end.lat, e2.end.lon);
1884 
1885  if ( check_intersection && edge_intersects(&A1, &A2, &B1, &B2) )
1886  {
1887  LWDEBUG(4,"edge intersection! returning 0.0");
1888  return 0.0;
1889  }
1890  d = s->radius * edge_distance_to_edge(&e1, &e2, &g1, &g2);
1891  LWDEBUGF(4,"got edge_distance_to_edge %.8g", d);
1892 
1893  if ( d < distance )
1894  {
1895  distance = d;
1896  nearest1 = g1;
1897  nearest2 = g2;
1898  }
1899  if ( d < tolerance )
1900  {
1901  if ( use_sphere )
1902  {
1903  return d;
1904  }
1905  else
1906  {
1907  d = spheroid_distance(&nearest1, &nearest2, s);
1908  if ( d < tolerance )
1909  return d;
1910  }
1911  }
1912 
1913  /* Copy end to start to allow a new end value in next iteration */
1914  e2.start = e2.end;
1915  B1 = B2;
1916  }
1917 
1918  /* Copy end to start to allow a new end value in next iteration */
1919  e1.start = e1.end;
1920  A1 = A2;
1921  LW_ON_INTERRUPT(return -1.0);
1922  }
1923  LWDEBUGF(4,"finished all loops, returning %.8g", distance);
1924 
1925  if ( use_sphere )
1926  return distance;
1927  else
1928  return spheroid_distance(&nearest1, &nearest2, s);
1929 }
1930 
1931 
1938 double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
1939 {
1940  int type;
1941  double radius2 = spheroid->radius * spheroid->radius;
1942 
1943  assert(lwgeom);
1944 
1945  /* No area in nothing */
1946  if ( lwgeom_is_empty(lwgeom) )
1947  return 0.0;
1948 
1949  /* Read the geometry type number */
1950  type = lwgeom->type;
1951 
1952  /* Anything but polygons and collections returns zero */
1953  if ( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
1954  return 0.0;
1955 
1956  /* Actually calculate area */
1957  if ( type == POLYGONTYPE )
1958  {
1959  LWPOLY *poly = (LWPOLY*)lwgeom;
1960  int i;
1961  double area = 0.0;
1962 
1963  /* Just in case there's no rings */
1964  if ( poly->nrings < 1 )
1965  return 0.0;
1966 
1967  /* First, the area of the outer ring */
1968  area += radius2 * ptarray_area_sphere(poly->rings[0]);
1969 
1970  /* Subtract areas of inner rings */
1971  for ( i = 1; i < poly->nrings; i++ )
1972  {
1973  area -= radius2 * ptarray_area_sphere(poly->rings[i]);
1974  }
1975  return area;
1976  }
1977 
1978  /* Recurse into sub-geometries to get area */
1979  if ( type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
1980  {
1981  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom;
1982  int i;
1983  double area = 0.0;
1984 
1985  for ( i = 0; i < col->ngeoms; i++ )
1986  {
1987  area += lwgeom_area_sphere(col->geoms[i], spheroid);
1988  }
1989  return area;
1990  }
1991 
1992  /* Shouldn't get here. */
1993  return 0.0;
1994 }
1995 
1996 
2006 LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, double distance, double azimuth)
2007 {
2008  GEOGRAPHIC_POINT geo_source, geo_dest;
2009  POINT4D pt_dest;
2010  double x, y;
2011  POINTARRAY *pa;
2012  LWPOINT *lwp;
2013 
2014  /* Check the azimuth validity, convert to radians */
2015  if ( azimuth < -2.0 * M_PI || azimuth > 2.0 * M_PI )
2016  {
2017  lwerror("Azimuth must be between -2PI and 2PI");
2018  return NULL;
2019  }
2020 
2021  /* Check the distance validity */
2022  if ( distance < 0.0 || distance > (M_PI * spheroid->radius) )
2023  {
2024  lwerror("Distance must be between 0 and %g", M_PI * spheroid->radius);
2025  return NULL;
2026  }
2027 
2028  /* Convert to ta geodetic point */
2029  x = lwpoint_get_x(r);
2030  y = lwpoint_get_y(r);
2031  geographic_point_init(x, y, &geo_source);
2032 
2033  /* Try the projection */
2034  if( spheroid_project(&geo_source, spheroid, distance, azimuth, &geo_dest) == LW_FAILURE )
2035  {
2036  LWDEBUGF(3, "Unable to project from (%g %g) with azimuth %g and distance %g", x, y, azimuth, distance);
2037  lwerror("Unable to project from (%g %g) with azimuth %g and distance %g", x, y, azimuth, distance);
2038  return NULL;
2039  }
2040 
2041  /* Build the output LWPOINT */
2042  pa = ptarray_construct(0, 0, 1);
2043  pt_dest.x = rad2deg(longitude_radians_normalize(geo_dest.lon));
2044  pt_dest.y = rad2deg(latitude_radians_normalize(geo_dest.lat));
2045  pt_dest.z = pt_dest.m = 0.0;
2046  ptarray_set_point4d(pa, 0, &pt_dest);
2047  lwp = lwpoint_construct(r->srid, NULL, pa);
2049  return lwp;
2050 }
2051 
2052 
2061 double lwgeom_azumith_spheroid(const LWPOINT *r, const LWPOINT *s, const SPHEROID *spheroid)
2062 {
2063  GEOGRAPHIC_POINT g1, g2;
2064  double x1, y1, x2, y2;
2065 
2066  /* Convert r to a geodetic point */
2067  x1 = lwpoint_get_x(r);
2068  y1 = lwpoint_get_y(r);
2069  geographic_point_init(x1, y1, &g1);
2070 
2071  /* Convert s to a geodetic point */
2072  x2 = lwpoint_get_x(s);
2073  y2 = lwpoint_get_y(s);
2074  geographic_point_init(x2, y2, &g2);
2075 
2076  /* Same point, return NaN */
2077  if ( FP_EQUALS(x1, x2) && FP_EQUALS(y1, y2) )
2078  {
2079  return NAN;
2080  }
2081 
2082  /* Do the direction calculation */
2083  return spheroid_direction(&g1, &g2, spheroid);
2084 }
2085 
2092 double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance)
2093 {
2094  uint8_t type1, type2;
2095  int check_intersection = LW_FALSE;
2096  GBOX gbox1, gbox2;
2097 
2098  gbox_init(&gbox1);
2099  gbox_init(&gbox2);
2100 
2101  assert(lwgeom1);
2102  assert(lwgeom2);
2103 
2104  LWDEBUGF(4, "entered function, tolerance %.8g", tolerance);
2105 
2106  /* What's the distance to an empty geometry? We don't know.
2107  Return a negative number so the caller can catch this case. */
2108  if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
2109  {
2110  return -1.0;
2111  }
2112 
2113  type1 = lwgeom1->type;
2114  type2 = lwgeom2->type;
2115 
2116  /* Make sure we have boxes */
2117  if ( lwgeom1->bbox )
2118  gbox1 = *(lwgeom1->bbox);
2119  else
2120  lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
2121 
2122  /* Make sure we have boxes */
2123  if ( lwgeom2->bbox )
2124  gbox2 = *(lwgeom2->bbox);
2125  else
2126  lwgeom_calculate_gbox_geodetic(lwgeom2, &gbox2);
2127 
2128  /* If the boxes aren't disjoint, we have to check for edge intersections */
2129  if ( gbox_overlaps(&gbox1, &gbox2) )
2130  check_intersection = LW_TRUE;
2131 
2132  /* Point/line combinations can all be handled with simple point array iterations */
2133  if ( ( type1 == POINTTYPE || type1 == LINETYPE ) &&
2134  ( type2 == POINTTYPE || type2 == LINETYPE ) )
2135  {
2136  POINTARRAY *pa1, *pa2;
2137 
2138  if ( type1 == POINTTYPE )
2139  pa1 = ((LWPOINT*)lwgeom1)->point;
2140  else
2141  pa1 = ((LWLINE*)lwgeom1)->points;
2142 
2143  if ( type2 == POINTTYPE )
2144  pa2 = ((LWPOINT*)lwgeom2)->point;
2145  else
2146  pa2 = ((LWLINE*)lwgeom2)->points;
2147 
2148  return ptarray_distance_spheroid(pa1, pa2, spheroid, tolerance, check_intersection);
2149  }
2150 
2151  /* Point/Polygon cases, if point-in-poly, return zero, else return distance. */
2152  if ( ( type1 == POLYGONTYPE && type2 == POINTTYPE ) ||
2153  ( type2 == POLYGONTYPE && type1 == POINTTYPE ) )
2154  {
2155  const POINT2D *p;
2156  LWPOLY *lwpoly;
2157  LWPOINT *lwpt;
2158  double distance = FLT_MAX;
2159  int i;
2160 
2161  if ( type1 == POINTTYPE )
2162  {
2163  lwpt = (LWPOINT*)lwgeom1;
2164  lwpoly = (LWPOLY*)lwgeom2;
2165  }
2166  else
2167  {
2168  lwpt = (LWPOINT*)lwgeom2;
2169  lwpoly = (LWPOLY*)lwgeom1;
2170  }
2171  p = getPoint2d_cp(lwpt->point, 0);
2172 
2173  /* Point in polygon implies zero distance */
2174  if ( lwpoly_covers_point2d(lwpoly, p) )
2175  {
2176  return 0.0;
2177  }
2178 
2179  /* Not inside, so what's the actual distance? */
2180  for ( i = 0; i < lwpoly->nrings; i++ )
2181  {
2182  double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwpt->point, spheroid, tolerance, check_intersection);
2183  if ( ring_distance < distance )
2184  distance = ring_distance;
2185  if ( distance < tolerance )
2186  return distance;
2187  }
2188  return distance;
2189  }
2190 
2191  /* Line/polygon case, if start point-in-poly, return zero, else return distance. */
2192  if ( ( type1 == POLYGONTYPE && type2 == LINETYPE ) ||
2193  ( type2 == POLYGONTYPE && type1 == LINETYPE ) )
2194  {
2195  const POINT2D *p;
2196  LWPOLY *lwpoly;
2197  LWLINE *lwline;
2198  double distance = FLT_MAX;
2199  int i;
2200 
2201  if ( type1 == LINETYPE )
2202  {
2203  lwline = (LWLINE*)lwgeom1;
2204  lwpoly = (LWPOLY*)lwgeom2;
2205  }
2206  else
2207  {
2208  lwline = (LWLINE*)lwgeom2;
2209  lwpoly = (LWPOLY*)lwgeom1;
2210  }
2211  p = getPoint2d_cp(lwline->points, 0);
2212 
2213  LWDEBUG(4, "checking if a point of line is in polygon");
2214 
2215  /* Point in polygon implies zero distance */
2216  if ( lwpoly_covers_point2d(lwpoly, p) )
2217  return 0.0;
2218 
2219  LWDEBUG(4, "checking ring distances");
2220 
2221  /* Not contained, so what's the actual distance? */
2222  for ( i = 0; i < lwpoly->nrings; i++ )
2223  {
2224  double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwline->points, spheroid, tolerance, check_intersection);
2225  LWDEBUGF(4, "ring[%d] ring_distance = %.8g", i, ring_distance);
2226  if ( ring_distance < distance )
2227  distance = ring_distance;
2228  if ( distance < tolerance )
2229  return distance;
2230  }
2231  LWDEBUGF(4, "all rings checked, returning distance = %.8g", distance);
2232  return distance;
2233 
2234  }
2235 
2236  /* Polygon/polygon case, if start point-in-poly, return zero, else return distance. */
2237  if ( ( type1 == POLYGONTYPE && type2 == POLYGONTYPE ) ||
2238  ( type2 == POLYGONTYPE && type1 == POLYGONTYPE ) )
2239  {
2240  const POINT2D *p;
2241  LWPOLY *lwpoly1 = (LWPOLY*)lwgeom1;
2242  LWPOLY *lwpoly2 = (LWPOLY*)lwgeom2;
2243  double distance = FLT_MAX;
2244  int i, j;
2245 
2246  /* Point of 2 in polygon 1 implies zero distance */
2247  p = getPoint2d_cp(lwpoly1->rings[0], 0);
2248  if ( lwpoly_covers_point2d(lwpoly2, p) )
2249  return 0.0;
2250 
2251  /* Point of 1 in polygon 2 implies zero distance */
2252  p = getPoint2d_cp(lwpoly2->rings[0], 0);
2253  if ( lwpoly_covers_point2d(lwpoly1, p) )
2254  return 0.0;
2255 
2256  /* Not contained, so what's the actual distance? */
2257  for ( i = 0; i < lwpoly1->nrings; i++ )
2258  {
2259  for ( j = 0; j < lwpoly2->nrings; j++ )
2260  {
2261  double ring_distance = ptarray_distance_spheroid(lwpoly1->rings[i], lwpoly2->rings[j], spheroid, tolerance, check_intersection);
2262  if ( ring_distance < distance )
2263  distance = ring_distance;
2264  if ( distance < tolerance )
2265  return distance;
2266  }
2267  }
2268  return distance;
2269  }
2270 
2271  /* Recurse into collections */
2272  if ( lwtype_is_collection(type1) )
2273  {
2274  int i;
2275  double distance = FLT_MAX;
2276  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1;
2277 
2278  for ( i = 0; i < col->ngeoms; i++ )
2279  {
2280  double geom_distance = lwgeom_distance_spheroid(col->geoms[i], lwgeom2, spheroid, tolerance);
2281  if ( geom_distance < distance )
2282  distance = geom_distance;
2283  if ( distance < tolerance )
2284  return distance;
2285  }
2286  return distance;
2287  }
2288 
2289  /* Recurse into collections */
2290  if ( lwtype_is_collection(type2) )
2291  {
2292  int i;
2293  double distance = FLT_MAX;
2294  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2;
2295 
2296  for ( i = 0; i < col->ngeoms; i++ )
2297  {
2298  double geom_distance = lwgeom_distance_spheroid(lwgeom1, col->geoms[i], spheroid, tolerance);
2299  if ( geom_distance < distance )
2300  distance = geom_distance;
2301  if ( distance < tolerance )
2302  return distance;
2303  }
2304  return distance;
2305  }
2306 
2307 
2308  lwerror("arguments include unsupported geometry type (%s, %s)", lwtype_name(type1), lwtype_name(type1));
2309  return -1.0;
2310 
2311 }
2312 
2313 
2314 int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
2315 {
2316  int type1, type2;
2317  GBOX gbox1, gbox2;
2318  gbox1.flags = gbox2.flags = 0;
2319 
2320  assert(lwgeom1);
2321  assert(lwgeom2);
2322 
2323  type1 = lwgeom1->type;
2324  type2 = lwgeom2->type;
2325 
2326  /* Currently a restricted implementation */
2327  if ( ! ( (type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE || type1 == COLLECTIONTYPE) &&
2328  (type2 == POINTTYPE || type2 == MULTIPOINTTYPE || type2 == COLLECTIONTYPE) ) )
2329  {
2330  lwerror("lwgeom_covers_lwgeom_sphere: only POLYGON covers POINT tests are currently supported");
2331  return LW_FALSE;
2332  }
2333 
2334  /* Make sure we have boxes */
2335  if ( lwgeom1->bbox )
2336  gbox1 = *(lwgeom1->bbox);
2337  else
2338  lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
2339 
2340  /* Make sure we have boxes */
2341  if ( lwgeom2->bbox )
2342  gbox2 = *(lwgeom2->bbox);
2343  else
2344  lwgeom_calculate_gbox_geodetic(lwgeom2, &gbox2);
2345 
2346 
2347  /* Handle the polygon/point case */
2348  if ( type1 == POLYGONTYPE && type2 == POINTTYPE )
2349  {
2350  POINT2D pt_to_test;
2351  getPoint2d_p(((LWPOINT*)lwgeom2)->point, 0, &pt_to_test);
2352  return lwpoly_covers_point2d((LWPOLY*)lwgeom1, &pt_to_test);
2353  }
2354 
2355  /* If any of the first argument parts covers the second argument, it's true */
2356  if ( lwtype_is_collection( type1 ) )
2357  {
2358  int i;
2359  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1;
2360 
2361  for ( i = 0; i < col->ngeoms; i++ )
2362  {
2363  if ( lwgeom_covers_lwgeom_sphere(col->geoms[i], lwgeom2) )
2364  {
2365  return LW_TRUE;
2366  }
2367  }
2368  return LW_FALSE;
2369  }
2370 
2371  /* Only if all of the second arguments are covered by the first argument is the condition true */
2372  if ( lwtype_is_collection( type2 ) )
2373  {
2374  int i;
2375  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2;
2376 
2377  for ( i = 0; i < col->ngeoms; i++ )
2378  {
2379  if ( ! lwgeom_covers_lwgeom_sphere(lwgeom1, col->geoms[i]) )
2380  {
2381  return LW_FALSE;
2382  }
2383  }
2384  return LW_TRUE;
2385  }
2386 
2387  /* Don't get here */
2388  lwerror("lwgeom_covers_lwgeom_sphere: reached end of function without resolution");
2389  return LW_FALSE;
2390 
2391 }
2392 
2398 int lwpoly_covers_point2d(const LWPOLY *poly, const POINT2D *pt_to_test)
2399 {
2400  int i;
2401  int in_hole_count = 0;
2402  POINT3D p;
2403  GEOGRAPHIC_POINT gpt_to_test;
2404  POINT2D pt_outside;
2405  GBOX gbox;
2406  gbox.flags = 0;
2407 
2408  /* Nulls and empties don't contain anything! */
2409  if ( ! poly || lwgeom_is_empty((LWGEOM*)poly) )
2410  {
2411  LWDEBUG(4,"returning false, geometry is empty or null");
2412  return LW_FALSE;
2413  }
2414 
2415  /* Make sure we have boxes */
2416  if ( poly->bbox )
2417  gbox = *(poly->bbox);
2418  else
2419  lwgeom_calculate_gbox_geodetic((LWGEOM*)poly, &gbox);
2420 
2421  /* Point not in box? Done! */
2422  geographic_point_init(pt_to_test->x, pt_to_test->y, &gpt_to_test);
2423  geog2cart(&gpt_to_test, &p);
2424  if ( ! gbox_contains_point3d(&gbox, &p) )
2425  {
2426  LWDEBUG(4, "the point is not in the box!");
2427  return LW_FALSE;
2428  }
2429 
2430  /* Calculate our outside point from the gbox */
2431  gbox_pt_outside(&gbox, &pt_outside);
2432 
2433  LWDEBUGF(4, "pt_outside POINT(%.18g %.18g)", pt_outside.x, pt_outside.y);
2434  LWDEBUGF(4, "pt_to_test POINT(%.18g %.18g)", pt_to_test->x, pt_to_test->y);
2435  LWDEBUGF(4, "polygon %s", lwgeom_to_ewkt((LWGEOM*)poly));
2436  LWDEBUGF(4, "gbox %s", gbox_to_string(&gbox));
2437 
2438  /* Not in outer ring? We're done! */
2439  if ( ! ptarray_contains_point_sphere(poly->rings[0], &pt_outside, pt_to_test) )
2440  {
2441  LWDEBUG(4,"returning false, point is outside ring");
2442  return LW_FALSE;
2443  }
2444 
2445  LWDEBUGF(4, "testing %d rings", poly->nrings);
2446 
2447  /* But maybe point is in a hole... */
2448  for ( i = 1; i < poly->nrings; i++ )
2449  {
2450  LWDEBUGF(4, "ring test loop %d", i);
2451  /* Count up hole containment. Odd => outside boundary. */
2452  if ( ptarray_contains_point_sphere(poly->rings[i], &pt_outside, pt_to_test) )
2453  in_hole_count++;
2454  }
2455 
2456  LWDEBUGF(4, "in_hole_count == %d", in_hole_count);
2457 
2458  if ( in_hole_count % 2 )
2459  {
2460  LWDEBUG(4,"returning false, inner ring containment count is odd");
2461  return LW_FALSE;
2462  }
2463 
2464  LWDEBUG(4,"returning true, inner ring containment count is even");
2465  return LW_TRUE;
2466 }
2467 
2468 
2473 int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point)
2474 {
2475  uint8_t *pa_ptr = NULL;
2476  assert(pa);
2477  assert(n >= 0);
2478  assert(n < pa->npoints);
2479 
2480  pa_ptr = getPoint_internal(pa, n);
2481  /* printf( "pa_ptr[0]: %g\n", *((double*)pa_ptr)); */
2482  *point = (POINT2D*)pa_ptr;
2483 
2484  return LW_SUCCESS;
2485 }
2486 
2488 {
2489  int i;
2490  int first = LW_TRUE;
2491  const POINT2D *p;
2492  POINT3D A1, A2;
2493  GBOX edge_gbox;
2494 
2495  assert(gbox);
2496  assert(pa);
2497 
2498  gbox_init(&edge_gbox);
2499  edge_gbox.flags = gbox->flags;
2500 
2501  if ( pa->npoints == 0 ) return LW_FAILURE;
2502 
2503  if ( pa->npoints == 1 )
2504  {
2505  p = getPoint2d_cp(pa, 0);
2506  ll2cart(p, &A1);
2507  gbox->xmin = gbox->xmax = A1.x;
2508  gbox->ymin = gbox->ymax = A1.y;
2509  gbox->zmin = gbox->zmax = A1.z;
2510  return LW_SUCCESS;
2511  }
2512 
2513  p = getPoint2d_cp(pa, 0);
2514  ll2cart(p, &A1);
2515 
2516  for ( i = 1; i < pa->npoints; i++ )
2517  {
2518 
2519  p = getPoint2d_cp(pa, i);
2520  ll2cart(p, &A2);
2521 
2522  edge_calculate_gbox(&A1, &A2, &edge_gbox);
2523 
2524  /* Initialize the box */
2525  if ( first )
2526  {
2527  gbox_duplicate(&edge_gbox, gbox);
2528  first = LW_FALSE;
2529  }
2530  /* Expand the box where necessary */
2531  else
2532  {
2533  gbox_merge(&edge_gbox, gbox);
2534  }
2535 
2536  A1 = A2;
2537  }
2538 
2539  return LW_SUCCESS;
2540 }
2541 
2542 static int lwpoint_calculate_gbox_geodetic(const LWPOINT *point, GBOX *gbox)
2543 {
2544  assert(point);
2545  return ptarray_calculate_gbox_geodetic(point->point, gbox);
2546 }
2547 
2548 static int lwline_calculate_gbox_geodetic(const LWLINE *line, GBOX *gbox)
2549 {
2550  assert(line);
2551  return ptarray_calculate_gbox_geodetic(line->points, gbox);
2552 }
2553 
2554 static int lwpolygon_calculate_gbox_geodetic(const LWPOLY *poly, GBOX *gbox)
2555 {
2556  GBOX ringbox;
2557  int i;
2558  int first = LW_TRUE;
2559  assert(poly);
2560  if ( poly->nrings == 0 )
2561  return LW_FAILURE;
2562  ringbox.flags = gbox->flags;
2563  for ( i = 0; i < poly->nrings; i++ )
2564  {
2565  if ( ptarray_calculate_gbox_geodetic(poly->rings[i], &ringbox) == LW_FAILURE )
2566  return LW_FAILURE;
2567  if ( first )
2568  {
2569  gbox_duplicate(&ringbox, gbox);
2570  first = LW_FALSE;
2571  }
2572  else
2573  {
2574  gbox_merge(&ringbox, gbox);
2575  }
2576  }
2577 
2578  /* If the box wraps a poly, push that axis to the absolute min/max as appropriate */
2579  gbox_check_poles(gbox);
2580 
2581  return LW_SUCCESS;
2582 }
2583 
2584 static int lwtriangle_calculate_gbox_geodetic(const LWTRIANGLE *triangle, GBOX *gbox)
2585 {
2586  assert(triangle);
2587  return ptarray_calculate_gbox_geodetic(triangle->points, gbox);
2588 }
2589 
2590 
2592 {
2593  GBOX subbox;
2594  int i;
2595  int result = LW_FAILURE;
2596  int first = LW_TRUE;
2597  assert(coll);
2598  if ( coll->ngeoms == 0 )
2599  return LW_FAILURE;
2600 
2601  subbox.flags = gbox->flags;
2602 
2603  for ( i = 0; i < coll->ngeoms; i++ )
2604  {
2605  if ( lwgeom_calculate_gbox_geodetic((LWGEOM*)(coll->geoms[i]), &subbox) == LW_SUCCESS )
2606  {
2607  /* Keep a copy of the sub-bounding box for later */
2608  if ( coll->geoms[i]->bbox )
2609  lwfree(coll->geoms[i]->bbox);
2610  coll->geoms[i]->bbox = gbox_copy(&subbox);
2611  if ( first )
2612  {
2613  gbox_duplicate(&subbox, gbox);
2614  first = LW_FALSE;
2615  }
2616  else
2617  {
2618  gbox_merge(&subbox, gbox);
2619  }
2620  result = LW_SUCCESS;
2621  }
2622  }
2623  return result;
2624 }
2625 
2627 {
2628  int result = LW_FAILURE;
2629  LWDEBUGF(4, "got type %d", geom->type);
2630 
2631  /* Add a geodetic flag to the incoming gbox */
2632  gbox->flags = gflags(FLAGS_GET_Z(geom->flags),FLAGS_GET_M(geom->flags),1);
2633 
2634  switch (geom->type)
2635  {
2636  case POINTTYPE:
2637  result = lwpoint_calculate_gbox_geodetic((LWPOINT*)geom, gbox);
2638  break;
2639  case LINETYPE:
2640  result = lwline_calculate_gbox_geodetic((LWLINE *)geom, gbox);
2641  break;
2642  case POLYGONTYPE:
2643  result = lwpolygon_calculate_gbox_geodetic((LWPOLY *)geom, gbox);
2644  break;
2645  case TRIANGLETYPE:
2646  result = lwtriangle_calculate_gbox_geodetic((LWTRIANGLE *)geom, gbox);
2647  break;
2648  case MULTIPOINTTYPE:
2649  case MULTILINETYPE:
2650  case MULTIPOLYGONTYPE:
2651  case POLYHEDRALSURFACETYPE:
2652  case TINTYPE:
2653  case COLLECTIONTYPE:
2654  result = lwcollection_calculate_gbox_geodetic((LWCOLLECTION *)geom, gbox);
2655  break;
2656  default:
2657  lwerror("lwgeom_calculate_gbox_geodetic: unsupported input geometry type: %d - %s",
2658  geom->type, lwtype_name(geom->type));
2659  break;
2660  }
2661  return result;
2662 }
2663 
2664 
2665 
2666 static int ptarray_check_geodetic(const POINTARRAY *pa)
2667 {
2668  int t;
2669  POINT2D pt;
2670 
2671  assert(pa);
2672 
2673  for (t=0; t<pa->npoints; t++)
2674  {
2675  getPoint2d_p(pa, t, &pt);
2676  /* printf( "%d (%g, %g)\n", t, pt.x, pt.y); */
2677  if ( pt.x < -180.0 || pt.y < -90.0 || pt.x > 180.0 || pt.y > 90.0 )
2678  return LW_FALSE;
2679  }
2680 
2681  return LW_TRUE;
2682 }
2683 
2684 static int lwpoint_check_geodetic(const LWPOINT *point)
2685 {
2686  assert(point);
2687  return ptarray_check_geodetic(point->point);
2688 }
2689 
2690 static int lwline_check_geodetic(const LWLINE *line)
2691 {
2692  assert(line);
2693  return ptarray_check_geodetic(line->points);
2694 }
2695 
2696 static int lwpoly_check_geodetic(const LWPOLY *poly)
2697 {
2698  int i = 0;
2699  assert(poly);
2700 
2701  for ( i = 0; i < poly->nrings; i++ )
2702  {
2703  if ( ptarray_check_geodetic(poly->rings[i]) == LW_FALSE )
2704  return LW_FALSE;
2705  }
2706  return LW_TRUE;
2707 }
2708 
2709 static int lwtriangle_check_geodetic(const LWTRIANGLE *triangle)
2710 {
2711  assert(triangle);
2712  return ptarray_check_geodetic(triangle->points);
2713 }
2714 
2715 
2717 {
2718  int i = 0;
2719  assert(col);
2720 
2721  for ( i = 0; i < col->ngeoms; i++ )
2722  {
2723  if ( lwgeom_check_geodetic(col->geoms[i]) == LW_FALSE )
2724  return LW_FALSE;
2725  }
2726  return LW_TRUE;
2727 }
2728 
2730 {
2731  if ( lwgeom_is_empty(geom) )
2732  return LW_TRUE;
2733 
2734  switch (geom->type)
2735  {
2736  case POINTTYPE:
2737  return lwpoint_check_geodetic((LWPOINT *)geom);
2738  case LINETYPE:
2739  return lwline_check_geodetic((LWLINE *)geom);
2740  case POLYGONTYPE:
2741  return lwpoly_check_geodetic((LWPOLY *)geom);
2742  case TRIANGLETYPE:
2743  return lwtriangle_check_geodetic((LWTRIANGLE *)geom);
2744  case MULTIPOINTTYPE:
2745  case MULTILINETYPE:
2746  case MULTIPOLYGONTYPE:
2747  case POLYHEDRALSURFACETYPE:
2748  case TINTYPE:
2749  case COLLECTIONTYPE:
2750  return lwcollection_check_geodetic((LWCOLLECTION *)geom);
2751  default:
2752  lwerror("lwgeom_check_geodetic: unsupported input geometry type: %d - %s",
2753  geom->type, lwtype_name(geom->type));
2754  }
2755  return LW_FALSE;
2756 }
2757 
2759 {
2760  int t;
2761  int changed = LW_FALSE;
2762  POINT4D pt;
2763 
2764  assert(pa);
2765 
2766  for ( t=0; t < pa->npoints; t++ )
2767  {
2768  getPoint4d_p(pa, t, &pt);
2769  if ( pt.x < -180.0 || pt.x > 180.0 || pt.y < -90.0 || pt.y > 90.0 )
2770  {
2771  pt.x = longitude_degrees_normalize(pt.x);
2772  pt.y = latitude_degrees_normalize(pt.y);
2773  ptarray_set_point4d(pa, t, &pt);
2774  changed = LW_TRUE;
2775  }
2776  }
2777  return changed;
2778 }
2779 
2781 {
2782  assert(point);
2783  return ptarray_force_geodetic(point->point);
2784 }
2785 
2787 {
2788  assert(line);
2789  return ptarray_force_geodetic(line->points);
2790 }
2791 
2793 {
2794  int i = 0;
2795  int changed = LW_FALSE;
2796  assert(poly);
2797 
2798  for ( i = 0; i < poly->nrings; i++ )
2799  {
2800  if ( ptarray_force_geodetic(poly->rings[i]) == LW_TRUE )
2801  changed = LW_TRUE;
2802  }
2803  return changed;
2804 }
2805 
2807 {
2808  int i = 0;
2809  int changed = LW_FALSE;
2810  assert(col);
2811 
2812  for ( i = 0; i < col->ngeoms; i++ )
2813  {
2814  if ( lwgeom_force_geodetic(col->geoms[i]) == LW_TRUE )
2815  changed = LW_TRUE;
2816  }
2817  return changed;
2818 }
2819 
2821 {
2822  switch ( lwgeom_get_type(geom) )
2823  {
2824  case POINTTYPE:
2825  return lwpoint_force_geodetic((LWPOINT *)geom);
2826  case LINETYPE:
2827  return lwline_force_geodetic((LWLINE *)geom);
2828  case POLYGONTYPE:
2829  return lwpoly_force_geodetic((LWPOLY *)geom);
2830  case MULTIPOINTTYPE:
2831  case MULTILINETYPE:
2832  case MULTIPOLYGONTYPE:
2833  case COLLECTIONTYPE:
2834  return lwcollection_force_geodetic((LWCOLLECTION *)geom);
2835  default:
2836  lwerror("unsupported input geometry type: %d", lwgeom_get_type(geom));
2837  }
2838  return LW_FALSE;
2839 }
2840 
2841 
2843 {
2844  GEOGRAPHIC_POINT a, b;
2845  double za = 0.0, zb = 0.0;
2846  POINT4D p;
2847  int i;
2848  int hasz = LW_FALSE;
2849  double length = 0.0;
2850  double seglength = 0.0;
2851 
2852  /* Return zero on non-sensical inputs */
2853  if ( ! pa || pa->npoints < 2 )
2854  return 0.0;
2855 
2856  /* See if we have a third dimension */
2857  hasz = FLAGS_GET_Z(pa->flags);
2858 
2859  /* Initialize first point */
2860  getPoint4d_p(pa, 0, &p);
2861  geographic_point_init(p.x, p.y, &a);
2862  if ( hasz )
2863  za = p.z;
2864 
2865  /* Loop and sum the length for each segment */
2866  for ( i = 1; i < pa->npoints; i++ )
2867  {
2868  seglength = 0.0;
2869  getPoint4d_p(pa, i, &p);
2870  geographic_point_init(p.x, p.y, &b);
2871  if ( hasz )
2872  zb = p.z;
2873 
2874  /* Special sphere case */
2875  if ( s->a == s->b )
2876  seglength = s->radius * sphere_distance(&a, &b);
2877  /* Spheroid case */
2878  else
2879  seglength = spheroid_distance(&a, &b, s);
2880 
2881  /* Add in the vertical displacement if we're in 3D */
2882  if ( hasz )
2883  seglength = sqrt( (zb-za)*(zb-za) + seglength*seglength );
2884 
2885  /* Add this segment length to the total */
2886  length += seglength;
2887 
2888  /* B gets incremented in the next loop, so we save the value here */
2889  a = b;
2890  za = zb;
2891  }
2892  return length;
2893 }
2894 
2895 double lwgeom_length_spheroid(const LWGEOM *geom, const SPHEROID *s)
2896 {
2897  int type;
2898  int i = 0;
2899  double length = 0.0;
2900 
2901  assert(geom);
2902 
2903  /* No area in nothing */
2904  if ( lwgeom_is_empty(geom) )
2905  return 0.0;
2906 
2907  type = geom->type;
2908 
2909  if ( type == POINTTYPE || type == MULTIPOINTTYPE )
2910  return 0.0;
2911 
2912  if ( type == LINETYPE )
2913  return ptarray_length_spheroid(((LWLINE*)geom)->points, s);
2914 
2915  if ( type == POLYGONTYPE )
2916  {
2917  LWPOLY *poly = (LWPOLY*)geom;
2918  for ( i = 0; i < poly->nrings; i++ )
2919  {
2920  length += ptarray_length_spheroid(poly->rings[i], s);
2921  }
2922  return length;
2923  }
2924 
2925  if ( type == TRIANGLETYPE )
2926  return ptarray_length_spheroid(((LWTRIANGLE*)geom)->points, s);
2927 
2928  if ( lwtype_is_collection( type ) )
2929  {
2930  LWCOLLECTION *col = (LWCOLLECTION*)geom;
2931 
2932  for ( i = 0; i < col->ngeoms; i++ )
2933  {
2934  length += lwgeom_length_spheroid(col->geoms[i], s);
2935  }
2936  return length;
2937  }
2938 
2939  lwerror("unsupported type passed to lwgeom_length_sphere");
2940  return 0.0;
2941 }
2942 
2949 static int
2951 {
2952 
2953  int i;
2954  POINT4D p;
2955  int altered = LW_FALSE;
2956  int rv = LW_FALSE;
2957  static double tolerance = 1e-10;
2958 
2959  if ( ! pa )
2960  lwerror("ptarray_nudge_geodetic called with null input");
2961 
2962  for(i = 0; i < pa->npoints; i++ )
2963  {
2964  getPoint4d_p(pa, i, &p);
2965  if ( p.x < -180.0 && (-180.0 - p.x < tolerance) )
2966  {
2967  p.x = -180.0;
2968  altered = LW_TRUE;
2969  }
2970  if ( p.x > 180.0 && (p.x - 180.0 < tolerance) )
2971  {
2972  p.x = 180.0;
2973  altered = LW_TRUE;
2974  }
2975  if ( p.y < -90.0 && (-90.0 - p.y < tolerance) )
2976  {
2977  p.y = -90.0;
2978  altered = LW_TRUE;
2979  }
2980  if ( p.y > 90.0 && (p.y - 90.0 < tolerance) )
2981  {
2982  p.y = 90.0;
2983  altered = LW_TRUE;
2984  }
2985  if ( altered == LW_TRUE )
2986  {
2987  ptarray_set_point4d(pa, i, &p);
2988  altered = LW_FALSE;
2989  rv = LW_TRUE;
2990  }
2991  }
2992  return rv;
2993 }
2994 
3001 int
3003 {
3004  int type;
3005  int i = 0;
3006  int rv = LW_FALSE;
3007 
3008  assert(geom);
3009 
3010  /* No points in nothing */
3011  if ( lwgeom_is_empty(geom) )
3012  return LW_FALSE;
3013 
3014  type = geom->type;
3015 
3016  if ( type == POINTTYPE )
3017  return ptarray_nudge_geodetic(((LWPOINT*)geom)->point);
3018 
3019  if ( type == LINETYPE )
3020  return ptarray_nudge_geodetic(((LWLINE*)geom)->points);
3021 
3022  if ( type == POLYGONTYPE )
3023  {
3024  LWPOLY *poly = (LWPOLY*)geom;
3025  for ( i = 0; i < poly->nrings; i++ )
3026  {
3027  int n = ptarray_nudge_geodetic(poly->rings[i]);
3028  rv = (rv == LW_TRUE ? rv : n);
3029  }
3030  return rv;
3031  }
3032 
3033  if ( type == TRIANGLETYPE )
3034  return ptarray_nudge_geodetic(((LWTRIANGLE*)geom)->points);
3035 
3036  if ( lwtype_is_collection( type ) )
3037  {
3038  LWCOLLECTION *col = (LWCOLLECTION*)geom;
3039 
3040  for ( i = 0; i < col->ngeoms; i++ )
3041  {
3042  int n = lwgeom_nudge_geodetic(col->geoms[i]);
3043  rv = (rv == LW_TRUE ? rv : n);
3044  }
3045  return rv;
3046  }
3047 
3048  lwerror("unsupported type (%s) passed to lwgeom_nudge_geodetic", lwtype_name(type));
3049  return rv;
3050 }
3051 
3052 
3056 static int
3057 point_in_cone(const POINT3D *A1, const POINT3D *A2, const POINT3D *P)
3058 {
3059  POINT3D AC; /* Center point of A1/A2 */
3060  double min_similarity, similarity;
3061 
3062  /* The normalized sum bisects the angle between start and end. */
3063  vector_sum(A1, A2, &AC);
3064  normalize(&AC);
3065 
3066  /* The projection of start onto the center defines the minimum similarity */
3067  min_similarity = dot_product(A1, &AC);
3068 
3069  /* The projection of candidate p onto the center */
3070  similarity = dot_product(P, &AC);
3071 
3072  /* If the point is more similar than the end, the point is in the cone */
3073  if ( similarity > min_similarity || fabs(similarity - min_similarity) < 2e-16 )
3074  {
3075  return LW_TRUE;
3076  }
3077  return LW_FALSE;
3078 }
3079 
3080 
3084 static int
3085 point3d_equals(const POINT3D *p1, const POINT3D *p2)
3086 {
3087  return FP_EQUALS(p1->x, p2->x) && FP_EQUALS(p1->y, p2->y) && FP_EQUALS(p1->z, p2->z);
3088 }
3089 
3094 static int
3095 dot_product_side(const POINT3D *p, const POINT3D *q)
3096 {
3097  double dp = dot_product(p, q);
3098 
3099  if ( FP_IS_ZERO(dp) )
3100  return 0;
3101 
3102  return dp < 0.0 ? -1 : 1;
3103 }
3104 
3109 int
3110 edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2)
3111 {
3112  POINT3D AN, BN, VN; /* Normals to plane A and plane B */
3113  double ab_dot;
3114  int a1_side, a2_side, b1_side, b2_side;
3115  int rv = PIR_NO_INTERACT;
3116 
3117  /* Normals to the A-plane and B-plane */
3118  unit_normal(A1, A2, &AN);
3119  unit_normal(B1, B2, &BN);
3120 
3121  /* Are A-plane and B-plane basically the same? */
3122  ab_dot = dot_product(&AN, &BN);
3123  if ( FP_EQUALS(fabs(ab_dot), 1.0) )
3124  {
3125  /* Co-linear case */
3126  if ( point_in_cone(A1, A2, B1) || point_in_cone(A1, A2, B2) ||
3127  point_in_cone(B1, B2, A1) || point_in_cone(B1, B2, A2) )
3128  {
3129  rv |= PIR_INTERSECTS;
3130  rv |= PIR_COLINEAR;
3131  }
3132  return rv;
3133  }
3134 
3135  /* What side of plane-A and plane-B do the end points */
3136  /* of A and B fall? */
3137  a1_side = dot_product_side(&BN, A1);
3138  a2_side = dot_product_side(&BN, A2);
3139  b1_side = dot_product_side(&AN, B1);
3140  b2_side = dot_product_side(&AN, B2);
3141 
3142  /* Both ends of A on the same side of plane B. */
3143  if ( a1_side == a2_side && a1_side != 0 )
3144  {
3145  /* No intersection. */
3146  return PIR_NO_INTERACT;
3147  }
3148 
3149  /* Both ends of B on the same side of plane A. */
3150  if ( b1_side == b2_side && b1_side != 0 )
3151  {
3152  /* No intersection. */
3153  return PIR_NO_INTERACT;
3154  }
3155 
3156  /* A straddles B and B straddles A, so... */
3157  if ( a1_side != a2_side && (a1_side + a2_side) == 0 &&
3158  b1_side != b2_side && (b1_side + b2_side) == 0 )
3159  {
3160  /* Have to check if intersection point is inside both arcs */
3161  unit_normal(&AN, &BN, &VN);
3162  if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
3163  {
3164  return PIR_INTERSECTS;
3165  }
3166 
3167  /* Have to check if intersection point is inside both arcs */
3168  vector_scale(&VN, -1);
3169  if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
3170  {
3171  return PIR_INTERSECTS;
3172  }
3173 
3174  return PIR_NO_INTERACT;
3175  }
3176 
3177  /* The rest are all intersects variants... */
3178  rv |= PIR_INTERSECTS;
3179 
3180  /* A touches B */
3181  if ( a1_side == 0 )
3182  {
3183  /* Touches at A1, A2 is on what side? */
3184  rv |= (a2_side < 0 ? PIR_A_TOUCH_RIGHT : PIR_A_TOUCH_LEFT);
3185  }
3186  else if ( a2_side == 0 )
3187  {
3188  /* Touches at A2, A1 is on what side? */
3189  rv |= (a1_side < 0 ? PIR_A_TOUCH_RIGHT : PIR_A_TOUCH_LEFT);
3190  }
3191 
3192  /* B touches A */
3193  if ( b1_side == 0 )
3194  {
3195  /* Touches at B1, B2 is on what side? */
3196  rv |= (b2_side < 0 ? PIR_B_TOUCH_RIGHT : PIR_B_TOUCH_LEFT);
3197  }
3198  else if ( b2_side == 0 )
3199  {
3200  /* Touches at B2, B1 is on what side? */
3201  rv |= (b1_side < 0 ? PIR_B_TOUCH_RIGHT : PIR_B_TOUCH_LEFT);
3202  }
3203 
3204  return rv;
3205 }
3206 
3215 int ptarray_contains_point_sphere(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test)
3216 {
3217  POINT3D S1, S2; /* Stab line end points */
3218  POINT3D E1, E2; /* Edge end points (3-space) */
3219  POINT2D p; /* Edge end points (lon/lat) */
3220  int count = 0, i, inter;
3221 
3222  /* Null input, not enough points for a ring? You ain't closed! */
3223  if ( ! pa || pa->npoints < 4 )
3224  return LW_FALSE;
3225 
3226  /* Set up our stab line */
3227  ll2cart(pt_to_test, &S1);
3228  ll2cart(pt_outside, &S2);
3229 
3230  /* Initialize first point */
3231  getPoint2d_p(pa, 0, &p);
3232  ll2cart(&p, &E1);
3233 
3234  /* Walk every edge and see if the stab line hits it */
3235  for ( i = 1; i < pa->npoints; i++ )
3236  {
3237  LWDEBUGF(4, "testing edge (%d)", i);
3238  LWDEBUGF(4, " start point == POINT(%.12g %.12g)", p.x, p.y);
3239 
3240  /* Read next point. */
3241  getPoint2d_p(pa, i, &p);
3242  ll2cart(&p, &E2);
3243 
3244  /* Skip over too-short edges. */
3245  if ( point3d_equals(&E1, &E2) )
3246  {
3247  continue;
3248  }
3249 
3250  /* Our test point is on an edge end! Point is "in ring" by our definition */
3251  if ( point3d_equals(&S1, &E1) )
3252  {
3253  return LW_TRUE;
3254  }
3255 
3256  /* Calculate relationship between stab line and edge */
3257  inter = edge_intersects(&S1, &S2, &E1, &E2);
3258 
3259  /* We have some kind of interaction... */
3260  if ( inter & PIR_INTERSECTS )
3261  {
3262  /* If the stabline is touching the edge, that implies the test point */
3263  /* is on the edge, so we're done, the point is in (on) the ring. */
3264  if ( (inter & PIR_A_TOUCH_RIGHT) || (inter & PIR_A_TOUCH_LEFT) )
3265  {
3266  return LW_TRUE;
3267  }
3268 
3269  /* It's a touching interaction, disregard all the left-side ones. */
3270  /* It's a co-linear intersection, ignore those. */
3271  if ( inter & PIR_B_TOUCH_RIGHT || inter & PIR_COLINEAR )
3272  {
3273  /* Do nothing, to avoid double counts. */
3274  LWDEBUGF(4," edge (%d) crossed, disregarding to avoid double count", i, count);
3275  }
3276  else
3277  {
3278  /* Increment crossingn count. */
3279  count++;
3280  LWDEBUGF(4," edge (%d) crossed, count == %d", i, count);
3281  }
3282  }
3283  else
3284  {
3285  LWDEBUGF(4," edge (%d) did not cross", i);
3286  }
3287 
3288  /* Increment to next edge */
3289  E1 = E2;
3290  }
3291 
3292  LWDEBUGF(4,"final count == %d", count);
3293 
3294  /* An odd number of crossings implies containment! */
3295  if ( count % 2 )
3296  {
3297  return LW_TRUE;
3298  }
3299 
3300  return LW_FALSE;
3301 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:549
double x
Definition: liblwgeom.h:351
#define LINETYPE
Definition: liblwgeom.h:85
GBOX * gbox_copy(const GBOX *box)
Return a copy of the GBOX, based on dimensionality of flags.
Definition: g_box.c:438
int gbox_centroid(const GBOX *gbox, POINT2D *out)
Computes the average(ish) center of the box and returns success.
Definition: lwgeodetic.c:258
void robust_cross_product(const GEOGRAPHIC_POINT *p, const GEOGRAPHIC_POINT *q, POINT3D *a)
Computes the cross product of two vectors using their lat, lng representations.
Definition: lwgeodetic.c:599
static int lwtriangle_calculate_gbox_geodetic(const LWTRIANGLE *triangle, GBOX *gbox)
Definition: lwgeodetic.c:2584
static int ptarray_nudge_geodetic(POINTARRAY *pa)
When features are snapped or sometimes they are just this way, they are very close to the geodetic bo...
Definition: lwgeodetic.c:2950
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
Definition: lwgeodetic.c:913
static int lwline_force_geodetic(LWLINE *line)
Definition: lwgeodetic.c:2786
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:62
double longitude_radians_normalize(double lon)
Convert a longitude to the range of -PI,PI.
Definition: lwgeodetic.c:41
GBOX * bbox
Definition: liblwgeom.h:397
POINTARRAY * points
Definition: liblwgeom.h:432
int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, GEOGRAPHIC_POINT *n)
Given a starting location r, a distance and an azimuth to the new point, compute the location of the ...
Definition: lwgeodetic.c:1279
double m
Definition: liblwgeom.h:351
double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2)
Calculate the distance between two edges.
Definition: lwgeodetic.c:1234
Two-point great circle segment from a to b.
Definition: lwgeodetic.h:56
double gbox_angular_height(const GBOX *gbox)
Returns the angular height (latitudinal span) of the box in radians.
Definition: lwgeodetic.c:179
#define PIR_B_TOUCH_RIGHT
Definition: lwgeodetic.h:86
double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid)
Computes the direction of the geodesic joining two points on the spheroid.
Definition: lwspheroid.c:282
void normalize(POINT3D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:580
static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY *pa2, const SPHEROID *s, double tolerance, int check_intersection)
Definition: lwgeodetic.c:1742
double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance)
Calculate the distance between two LWGEOMs, using the coordinates are longitude and latitude...
Definition: lwgeodetic.c:2092
char * r
Definition: cu_in_wkt.c:24
int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point)
This function can only be used on LWGEOM that is built on top of GSERIALIZED, otherwise alignment err...
Definition: lwgeodetic.c:2473
void gbox_duplicate(const GBOX *original, GBOX *duplicate)
Copy the values of original GBOX into duplicate.
Definition: g_box.c:445
void lwfree(void *mem)
Definition: lwutil.c:242
int edge_contains_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
Returns true if the point p is on the minor edge defined by the end points of e.
Definition: lwgeodetic.c:999
double y
Definition: liblwgeom.h:339
char * gbox_to_string(const GBOX *gbox)
Allocate a string representation of the GBOX, based on dimensionality of flags.
Definition: g_box.c:404
int npoints
Definition: liblwgeom.h:370
static int point3d_equals(const POINT3D *p1, const POINT3D *p2)
Utility function for ptarray_contains_point_sphere()
Definition: lwgeodetic.c:3085
uint8_t type
Definition: liblwgeom.h:502
int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Definition: lwgeodetic.c:631
#define PIR_A_TOUCH_RIGHT
Definition: lwgeodetic.h:84
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:842
int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
Calculate covers predicate for two lwgeoms on the sphere.
Definition: lwgeodetic.c:2314
double z_to_latitude(double z, int top)
Used in great circle to compute the pole of the great circle.
Definition: lwgeodetic.c:1014
#define POLYGONTYPE
Definition: liblwgeom.h:86
Datum area(PG_FUNCTION_ARGS)
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:482
double b
Definition: liblwgeom.h:313
uint8_t flags
Definition: liblwgeom.h:396
double xmax
Definition: liblwgeom.h:292
#define NAN
Definition: lwgeodetic.h:37
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
int gbox_contains_point3d(const GBOX *gbox, const POINT3D *pt)
Return true if the point is inside the gbox.
Definition: g_box.c:259
#define LW_SUCCESS
Definition: liblwgeom.h:79
void cart2geog(const POINT3D *p, GEOGRAPHIC_POINT *g)
Convert cartesion coordinates on unit sphere to spherical coordinates.
Definition: lwgeodetic.c:379
char * w
Definition: cu_out_twkb.c:25
double radius
Definition: liblwgeom.h:317
int edge_calculate_gbox_slow(const GEOGRAPHIC_EDGE *e, GBOX *gbox)
Definition: lwgeodetic.c:1309
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define TRIANGLETYPE
Definition: liblwgeom.h:97
static int dot_product_side(const POINT3D *p, const POINT3D *q)
Utility function for edge_intersects(), signum with a tolerance in determining if the value is zero...
Definition: lwgeodetic.c:3095
#define LW_ON_INTERRUPT(x)
double x
Definition: liblwgeom.h:339
#define POLYHEDRALSURFACETYPE
Definition: liblwgeom.h:96
void point_shift(GEOGRAPHIC_POINT *p, double shift)
Shift a point around by a number of radians.
Definition: lwgeodetic.c:151
void y_to_z(POINT3D *p)
Definition: lwgeodetic.c:623
double lwgeom_length_spheroid(const LWGEOM *geom, const SPHEROID *s)
Calculate the geodetic length of a lwgeom on the unit sphere.
Definition: lwgeodetic.c:2895
static void vector_scale(POINT3D *n, double scale)
Scale a vector out by a factor.
Definition: lwgeodetic.c:452
int gbox_init_point3d(const POINT3D *p, GBOX *gbox)
Initialize a GBOX using the values of the point.
Definition: g_box.c:251
int lwgeom_force_geodetic(LWGEOM *geom)
Force coordinates of LWGEOM into geodetic range (-180, -90, 180, 90)
Definition: lwgeodetic.c:2820
static int lwline_calculate_gbox_geodetic(const LWLINE *line, GBOX *gbox)
Definition: lwgeodetic.c:2548
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:446
GBOX * bbox
Definition: liblwgeom.h:452
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:138
double longitude_degrees_normalize(double lon)
Convert a longitude to the range of -180,180.
Definition: lwgeodetic.c:97
#define FP_MIN(A, B)
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:47
POINTARRAY * point
Definition: liblwgeom.h:410
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:849
int32_t srid
Definition: liblwgeom.h:398
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:252
#define FP_IS_ZERO(A)
static double sphere_angle(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c)
Returns the angle in radians at point B of the triangle formed by A-B-C.
Definition: lwgeodetic.c:686
int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
Given a location, an azimuth and a distance, computes the location of the projected point...
Definition: lwspheroid.c:359
int gbox_merge(const GBOX *new_box, GBOX *merge_box)
Update the merged GBOX to be large enough to include itself and the new box.
Definition: g_box.c:269
#define PIR_A_TOUCH_LEFT
Definition: lwgeodetic.h:85
#define LW_FAILURE
Definition: liblwgeom.h:78
double latitude_degrees_normalize(double lat)
Convert a latitude to the range of -90,90.
Definition: lwgeodetic.c:124
double z
Definition: liblwgeom.h:339
double x
Definition: liblwgeom.h:327
int clairaut_geographic(const GEOGRAPHIC_POINT *start, const GEOGRAPHIC_POINT *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
Computes the pole of the great circle disk which is the intersection of the great circle with the lin...
Definition: lwgeodetic.c:1066
void unit_normal(const POINT3D *P1, const POINT3D *P2, POINT3D *normal)
Calculates the unit normal to two vectors, trying to avoid problems with over-narrow or over-wide cas...
Definition: lwgeodetic.c:506
double lwpoint_get_x(const LWPOINT *point)
Definition: lwpoint.c:63
int ptarray_contains_point_sphere(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test)
This routine returns LW_TRUE if the stabline joining the pt_outside and pt_to_test crosses the ring a...
Definition: lwgeodetic.c:3215
static int lwpoint_force_geodetic(LWPOINT *point)
Definition: lwgeodetic.c:2780
double zmax
Definition: liblwgeom.h:296
double ymin
Definition: liblwgeom.h:293
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:216
static double dot_product(const POINT3D *p1, const POINT3D *p2)
Convert cartesion coordinates on unit sphere to lon/lat coordinates static void cart2ll(const POINT3D...
Definition: lwgeodetic.c:411
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:262
void gbox_init(GBOX *gbox)
Zero out all the entries in the GBOX.
Definition: g_box.c:51
double xmin
Definition: liblwgeom.h:291
int ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int allow_duplicates)
Append a point to the end of an existing POINTARRAY If allow_duplicate is LW_FALSE, then a duplicate point will not be added.
Definition: ptarray.c:156
#define LW_FALSE
Definition: liblwgeom.h:76
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, int n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:485
uint8_t flags
Definition: liblwgeom.h:368
#define rad2deg(r)
Definition: lwgeodetic.h:75
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:58
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
Definition: lwgeodetic.c:1183
int edge_point_in_cone(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
Returns true if the point p is inside the cone defined by the two ends of the edge e...
Definition: lwgeodetic.c:753
double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the area of an LWGEOM.
Definition: lwgeodetic.c:1938
LWGEOM ** geoms
Definition: liblwgeom.h:508
double ptarray_length_spheroid(const POINTARRAY *pa, const SPHEROID *s)
Definition: lwgeodetic.c:2842
int edge_contains_coplanar_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
True if the longitude of p is within the range of the longitude of the ends of e. ...
Definition: lwgeodetic.c:800
void ll2cart(const POINT2D *g, POINT3D *p)
Convert lon/lat coordinates to cartesion coordinates on unit sphere.
Definition: lwgeodetic.c:388
int clairaut_cartesian(const POINT3D *start, const POINT3D *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
Computes the pole of the great circle disk which is the intersection of the great circle with the lin...
Definition: lwgeodetic.c:1041
#define TINTYPE
Definition: liblwgeom.h:98
void vector_rotate(const POINT3D *v1, const POINT3D *v2, double angle, POINT3D *n)
Rotates v1 through an angle (in radians) within the plane defined by v1/v2, returns the rotated vecto...
Definition: lwgeodetic.c:538
void x_to_z(POINT3D *p)
Definition: lwgeodetic.c:616
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:2626
int gbox_merge_point3d(const POINT3D *p, GBOX *gbox)
Update the GBOX to be large enough to include itself and the new point.
Definition: g_box.c:240
static void vector_difference(const POINT3D *a, const POINT3D *b, POINT3D *n)
Calculate the difference of two vectors.
Definition: lwgeodetic.c:441
uint8_t * getPoint_internal(const POINTARRAY *pa, int n)
Definition: ptarray.c:1713
static int lwpoint_calculate_gbox_geodetic(const LWPOINT *point, GBOX *gbox)
Definition: lwgeodetic.c:2542
static int lwcollection_force_geodetic(LWCOLLECTION *col)
Definition: lwgeodetic.c:2806
POINTARRAY ** rings
Definition: liblwgeom.h:456
int count
Definition: genraster.py:56
int lwtype_is_collection(uint8_t type)
Determine whether a type number is a collection or not.
Definition: lwgeom.c:1012
static POINTARRAY * ptarray_segmentize_sphere(const POINTARRAY *pa_in, double max_seg_length)
Create a new point array with no segment longer than the input segment length (expressed in radians!)...
Definition: lwgeodetic.c:1591
void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside)
Definition: lwgeodetic.c:1440
int nrings
Definition: liblwgeom.h:454
double ymax
Definition: liblwgeom.h:294
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:59
int lwgeom_nudge_geodetic(LWGEOM *geom)
When features are snapped or sometimes they are just this way, they are very close to the geodetic bo...
Definition: lwgeodetic.c:3002
char * s
Definition: cu_in_wkt.c:23
double y
Definition: liblwgeom.h:327
static int point_in_cone(const POINT3D *A1, const POINT3D *A2, const POINT3D *P)
Utility function for checking if P is within the cone defined by A1/A2.
Definition: lwgeodetic.c:3057
int ptarray_calculate_gbox_geodetic(const POINTARRAY *pa, GBOX *gbox)
Calculate geodetic (x/y/z) box and add values to gbox.
Definition: lwgeodetic.c:2487
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:461
#define FLAGS_GET_Z(flags)
Macros for manipulating the &#39;flags&#39; byte.
Definition: liblwgeom.h:139
double gbox_angular_width(const GBOX *gbox)
Returns the angular width (longitudinal span) of the box in radians.
Definition: lwgeodetic.c:206
double z
Definition: liblwgeom.h:351
static int lwpolygon_calculate_gbox_geodetic(const LWPOLY *poly, GBOX *gbox)
Definition: lwgeodetic.c:2554
#define POW2(x)
Definition: lwgeodetic.h:42
LWGEOM * lwgeom_clone(const LWGEOM *lwgeom)
Clone LWGEOM object.
Definition: lwgeom.c:408
Datum distance(PG_FUNCTION_ARGS)
static int lwline_check_geodetic(const LWLINE *line)
Definition: lwgeodetic.c:2690
void lwgeom_set_geodetic(LWGEOM *geom, int value)
Set the FLAGS geodetic bit on geometry an all sub-geometries and pointlists.
Definition: lwgeom.c:871
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:102
int32_t srid
Definition: liblwgeom.h:409
uint8_t flags
Definition: liblwgeom.h:290
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:145
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesion coordinates on unit sphere.
Definition: lwgeodetic.c:369
void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
Given a unit geocentric gbox, return a lon/lat (degrees) coordinate point point that is guaranteed to...
Definition: lwgeodetic.c:1461
int gbox_geocentric_slow
For testing geodetic bounding box, we have a magic global variable.
Definition: lwgeodetic.c:36
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:89
void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n)
Calculate the sum of two vectors.
Definition: lwgeodetic.c:430
double latitude_radians_normalize(double lat)
Convert a latitude to the range of -PI/2,PI/2.
Definition: lwgeodetic.c:69
int edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2)
Returns non-zero if edges A and B interact.
Definition: lwgeodetic.c:3110
double sphere_distance_cartesian(const POINT3D *s, const POINT3D *e)
Given two unit vectors, calculate their distance apart in radians.
Definition: lwgeodetic.c:932
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:43
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:161
LWGEOM * lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length)
Create a new, densified geometry where no segment is longer than max_seg_length.
Definition: lwgeodetic.c:1650
double a
Definition: liblwgeom.h:312
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:156
static int lwcollection_check_geodetic(const LWCOLLECTION *col)
Definition: lwgeodetic.c:2716
double zmin
Definition: liblwgeom.h:295
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:171
double ptarray_area_sphere(const POINTARRAY *pa)
Returns the area of the ring (ring must be closed) in square radians (surface of the sphere is 4*PI)...
Definition: lwgeodetic.c:1714
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
static int lwtriangle_check_geodetic(const LWTRIANGLE *triangle)
Definition: lwgeodetic.c:2709
double lwgeom_azumith_spheroid(const LWPOINT *r, const LWPOINT *s, const SPHEROID *spheroid)
Calculate a bearing (azimuth) given a source and destination point.
Definition: lwgeodetic.c:2061
LWPOLY * lwpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoly.c:139
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:140
int p3d_same(const POINT3D *p1, const POINT3D *p2)
Definition: lwalgorithm.c:40
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:74
int gbox_overlaps(const GBOX *g1, const GBOX *g2)
Return LW_TRUE if the GBOX overlaps, LW_FALSE otherwise.
Definition: g_box.c:295
uint8_t type
Definition: liblwgeom.h:395
static double sphere_signed_area(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c)
Computes the spherical area of a triangle.
Definition: lwgeodetic.c:706
type
Definition: ovdump.py:41
#define FP_EQUALS(A, B)
double lwpoint_get_y(const LWPOINT *point)
Definition: lwpoint.c:73
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:267
LWPOINT * lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, double distance, double azimuth)
Calculate a projected point given a source point, a distance and a bearing.
Definition: lwgeodetic.c:2006
static int lwpoly_check_geodetic(const LWPOLY *poly)
Definition: lwgeodetic.c:2696
double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
Computes the shortest distance along the surface of the spheroid between two points.
Definition: lwspheroid.c:186
int lwpoly_add_ring(LWPOLY *poly, POINTARRAY *pa)
Add a ring, allocating extra space if necessary.
Definition: lwpoly.c:227
static int lwpoly_force_geodetic(LWPOLY *poly)
Definition: lwgeodetic.c:2792
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:36
#define PIR_NO_INTERACT
Bitmask elements for edge_intersects() return value.
Definition: lwgeodetic.h:81
int lwpoly_covers_point2d(const LWPOLY *poly, const POINT2D *pt_to_test)
Given a polygon (lon/lat decimal degrees) and point (lon/lat decimal degrees) and a guaranteed outsid...
Definition: lwgeodetic.c:2398
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:111
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1310
#define PIR_B_TOUCH_LEFT
Definition: lwgeodetic.h:87
static int lwpoint_check_geodetic(const LWPOINT *point)
Definition: lwgeodetic.c:2684
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:64
double y
Definition: liblwgeom.h:351
#define MULTILINETYPE
Definition: liblwgeom.h:88
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
Returns true if the point p is on the great circle plane.
Definition: lwgeodetic.c:740
double vector_angle(const POINT3D *v1, const POINT3D *v2)
Angle between two unit vectors.
Definition: lwgeodetic.c:470
static int gbox_check_poles(GBOX *gbox)
Check to see if this geocentric gbox is wrapped around a pole.
Definition: lwgeodetic.c:307
static int edge_point_side(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
Returns -1 if the point is to the left of the plane formed by the edge, 1 if the point is to the righ...
Definition: lwgeodetic.c:659
static void cross_product(const POINT3D *a, const POINT3D *b, POINT3D *n)
Calculate the cross product of two vectors.
Definition: lwgeodetic.c:419
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
static void normalize2d(POINT2D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:489
static int ptarray_force_geodetic(POINTARRAY *pa)
Definition: lwgeodetic.c:2758
#define PIR_COLINEAR
Definition: lwgeodetic.h:83
int edge_calculate_gbox(const POINT3D *A1, const POINT3D *A2, GBOX *gbox)
The magic function, given an edge in spherical coordinates, calculate a 3D bounding box that fully co...
Definition: lwgeodetic.c:1373
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:856
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:187
static int ptarray_segmentize_sphere_edge_recursive(const POINT3D *p1, const POINT3D *p2, const POINT4D *v1, const POINT4D *v2, double d, double max_seg_length, POINTARRAY *pa)
Definition: lwgeodetic.c:1542
int lwgeom_check_geodetic(const LWGEOM *geom)
Check that coordinates of LWGEOM are all within the geodetic range (-180, -90, 180, 90)
Definition: lwgeodetic.c:2729
int edge_intersection(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *g)
Returns true if an intersection can be calculated, and places it in *g.
Definition: lwgeodetic.c:1092
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:102
static int lwcollection_calculate_gbox_geodetic(const LWCOLLECTION *coll, GBOX *gbox)
Definition: lwgeodetic.c:2591
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:244
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
int p4d_same(const POINT4D *p1, const POINT4D *p2)
Definition: lwalgorithm.c:31
#define SIGNUM(n)
Macro that returns: -1 if n < 0, 1 if n > 0, 0 if n == 0.
#define FP_MAX(A, B)
#define PIR_INTERSECTS
Definition: lwgeodetic.h:82
POINTARRAY * points
Definition: liblwgeom.h:421
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:232
double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, double d)
Given two points on a unit sphere, calculate the direction from s to e.
Definition: lwgeodetic.c:940
static int ptarray_check_geodetic(const POINTARRAY *pa)
Definition: lwgeodetic.c:2666