PostGIS  2.1.10dev-r@@SVN_REVISION@@
cu_geodetic.c
Go to the documentation of this file.
1 /**********************************************************************
2  * $Id: cu_geodetic.c 14487 2015-12-14 13:10:24Z dbaston $
3  *
4  * PostGIS - Spatial Types for PostgreSQL
5  * http://postgis.net
6  * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <string.h>
16 #include "CUnit/Basic.h"
17 
18 #include "liblwgeom_internal.h"
19 #include "lwgeodetic.h"
20 #include "cu_tester.h"
21 
22 #define RANDOM_TEST 0
23 
28 {
29  (e->start).lat = deg2rad((e->start).lat);
30  (e->end).lat = deg2rad((e->end).lat);
31  (e->start).lon = deg2rad((e->start).lon);
32  (e->end).lon = deg2rad((e->end).lon);
33 }
34 
50 {
53 }
54 
59 {
60  p->lat = rad2deg(p->lat);
61  p->lon = rad2deg(p->lon);
62 }
63 
64 static void test_signum(void)
65 {
66  CU_ASSERT_EQUAL(signum(-5.0),-1);
67  CU_ASSERT_EQUAL(signum(5.0),1);
68 }
69 
70 
71 static void test_sphere_direction(void)
72 {
74  double dir, dist;
75 
76  geographic_point_init(0, 0, &s);
77  geographic_point_init(1, 0, &e);
78  dist = sphere_distance(&s, &e);
79  dir = sphere_direction(&s, &e, dist);
80  CU_ASSERT_DOUBLE_EQUAL(dir, M_PI / 2.0, 0.0001);
81 
82  geographic_point_init(0, 0, &s);
83  geographic_point_init(0, 1, &e);
84  dist = sphere_distance(&s, &e);
85  dir = sphere_direction(&s, &e, dist);
86  CU_ASSERT_DOUBLE_EQUAL(dir, 0.0, 0.0001);
87 
88 }
89 
90 static void test_sphere_project(void)
91 {
93  double dir1, dist1, dir2, dist2;
94 
95  dir1 = M_PI/2;
96  dist1 = 0.1;
97 
98  geographic_point_init(0, 0, &s);
99  sphere_project(&s, dist1, dir1, &e);
100 
101  dist2 = sphere_distance(&s, &e);
102  dir2 = sphere_direction(&s, &e, dist1);
103 
104  CU_ASSERT_DOUBLE_EQUAL(dist1, dist2, 0.0001);
105  CU_ASSERT_DOUBLE_EQUAL(dir1, dir2, 0.0001);
106 
107  dist1 = sphere_distance(&e, &s);
108  dir1 = sphere_direction(&e, &s, dist1);
109  sphere_project(&e, dist1, dir1, &s);
110 
111  CU_ASSERT_DOUBLE_EQUAL(s.lon, 0.0, 0.0001);
112  CU_ASSERT_DOUBLE_EQUAL(s.lat, 0.0, 0.0001);
113 
114  geographic_point_init(0, 0.2, &e);
115  geographic_point_init(0, 0.4, &s);
116  dist1 = sphere_distance(&s, &e);
117  dir1 = sphere_direction(&e, &s, dist1);
118  CU_ASSERT_DOUBLE_EQUAL(dir1, 0.0, 0.0001);
119 
120  geographic_point_init(0, 1, &s);
121  geographic_point_init(0, 2, &e);
122  dist2 = sphere_distance(&s, &e);
123  dir2 = sphere_direction(&s, &e, dist2);
124  CU_ASSERT_DOUBLE_EQUAL(dir2, 0.0, 0.0001);
125 
126  geographic_point_init(1, 1, &e);
127  dist2 = sphere_distance(&s, &e);
128  dir2 = sphere_direction(&s, &e, dist2);
129  CU_ASSERT_DOUBLE_EQUAL(dir2, 1.57064, 0.0001);
130 
131  geographic_point_init(0, 0, &e);
132  dist2 = sphere_distance(&s, &e);
133  dir2 = sphere_direction(&s, &e, dist2);
134  CU_ASSERT_DOUBLE_EQUAL(dir2, 3.14159, 0.0001);
135 
136  geographic_point_init(-1, 1, &e);
137  dist2 = sphere_distance(&s, &e);
138  dir2 = sphere_direction(&s, &e, dist2);
139  CU_ASSERT_DOUBLE_EQUAL(dir2, -1.57064, 0.0001);
140 
141  geographic_point_init(1, 2, &e);
142  dist2 = sphere_distance(&s, &e);
143  dir2 = sphere_direction(&s, &e, dist2);
144  CU_ASSERT_DOUBLE_EQUAL(dir2, 0.785017, 0.0001);
145 
146  geographic_point_init(-1, 0, &e);
147  dist2 = sphere_distance(&s, &e);
148  dir2 = sphere_direction(&s, &e, dist2);
149  CU_ASSERT_DOUBLE_EQUAL(dir2, -2.35612, 0.0001);
150 }
151 
152 #if 0
153 
157 static void cross_product_stability(void)
158 {
159  POINT2D p1, p2;
160  int i;
161  GEOGRAPHIC_POINT g1, g2;
162  POINT3D A1, A2;
163  POINT3D Nr, Nc;
164  POINT3D Or, Oc;
165 
166  p1.x = 10.0;
167  p1.y = 45.0;
168  p2.x = 10.0;
169  p2.y = 50.0;
170 
171  geographic_point_init(p1.x, p1.y, &g1);
172  ll2cart(&p1, &A1);
173 
174  for ( i = 0; i < 40; i++ )
175  {
176  geographic_point_init(p2.x, p2.y, &g2);
177  ll2cart(&p2, &A2);
178 
179  /* Skea */
180  robust_cross_product(&g1, &g2, &Nr);
181  normalize(&Nr);
182 
183  /* Ramsey */
184  unit_normal(&A1, &A2, &Nc);
185 
186  if ( i > 0 )
187  {
188  printf("\n- %d -------------------- %.24g ------------------------\n", i, p2.y);
189  printf("Skea: %.24g,%.24g,%.24g\n", Nr.x, Nr.y, Nr.z);
190  printf("Skea Diff: %.24g,%.24g,%.24g\n", Or.x-Nr.x, Or.y-Nr.y, Or.z-Nr.z);
191  printf("Ramsey: %.24g,%.24g,%.24g\n", Nc.x, Nc.y, Nc.z);
192  printf("Ramsey Diff: %.24g,%.24g,%.24g\n", Oc.x-Nc.x, Oc.y-Nc.y, Oc.z-Nc.z);
193  printf("Diff: %.24g,%.24g,%.24g\n", Nr.x-Nc.x, Nr.y-Nc.y, Nr.z-Nc.z);
194  }
195 
196  Or = Nr;
197  Oc = Nc;
198 
199  p2.y += (p1.y - p2.y)/2.0;
200  }
201 }
202 #endif
203 
205 {
206 #if RANDOM_TEST
207  const double gtolerance = 0.000001;
208  const int loops = RANDOM_TEST;
209  int i;
210  double ll[64];
211  GBOX gbox;
212  GBOX gbox_slow;
213  int rndlat;
214  int rndlon;
215 
216  POINTARRAY *pa;
217  LWGEOM *lwline;
218 
219  ll[0] = -3.083333333333333333333333333333333;
220  ll[1] = 9.83333333333333333333333333333333;
221  ll[2] = 15.5;
222  ll[3] = -5.25;
223 
224  pa = ptarray_construct_reference_data(0, 0, 2, (uint8_t*)ll);
225 
227  FLAGS_SET_GEODETIC(lwline->flags, 1);
228 
229  srandomdev();
230 
231  for ( i = 0; i < loops; i++ )
232  {
233  rndlat = (int)(90.0 - 180.0 * (double)random() / pow(2.0, 31.0));
234  rndlon = (int)(180.0 - 360.0 * (double)random() / pow(2.0, 31.0));
235  ll[0] = (double)rndlon;
236  ll[1] = (double)rndlat;
237 
238  rndlat = (int)(90.0 - 180.0 * (double)random() / pow(2.0, 31.0));
239  rndlon = (int)(180.0 - 360.0 * (double)random() / pow(2.0, 31.0));
240  ll[2] = (double)rndlon;
241  ll[3] = (double)rndlat;
242 
244  lwgeom_calculate_gbox_geodetic(lwline, &gbox);
246  lwgeom_calculate_gbox_geodetic(lwline, &gbox_slow);
248 
249  if (
250  ( fabs( gbox.xmin - gbox_slow.xmin ) > gtolerance ) ||
251  ( fabs( gbox.xmax - gbox_slow.xmax ) > gtolerance ) ||
252  ( fabs( gbox.ymin - gbox_slow.ymin ) > gtolerance ) ||
253  ( fabs( gbox.ymax - gbox_slow.ymax ) > gtolerance ) ||
254  ( fabs( gbox.zmin - gbox_slow.zmin ) > gtolerance ) ||
255  ( fabs( gbox.zmax - gbox_slow.zmax ) > gtolerance ) )
256  {
257  printf("\n-------\n");
258  printf("If you are seeing this, cut and paste it, it is a randomly generated test case!\n");
259  printf("LOOP: %d\n", i);
260  printf("SEGMENT (Lon Lat): (%.9g %.9g) (%.9g %.9g)\n", ll[0], ll[1], ll[2], ll[3]);
261  printf("CALC: %s\n", gbox_to_string(&gbox));
262  printf("SLOW: %s\n", gbox_to_string(&gbox_slow));
263  printf("-------\n\n");
264  CU_FAIL_FATAL(Slow (GOOD) and fast (CALC) box calculations returned different values!!);
265  }
266 
267  }
268 
269  lwgeom_free(lwline);
270 #endif /* RANDOM_TEST */
271 }
272 
273 #include "cu_geodetic_data.h"
274 
276 {
277  LWGEOM *lwg;
278  GBOX gbox, gbox_slow;
279  int i;
280 
281  for ( i = 0; i < gbox_data_length; i++ )
282  {
283 #if 0
284 // if ( i != 0 ) continue; /* skip our bad case */
285  printf("\n\n------------\n");
286  printf("%s\n", gbox_data[i]);
287 #endif
289  FLAGS_SET_GEODETIC(lwg->flags, 1);
291  lwgeom_calculate_gbox(lwg, &gbox);
293  lwgeom_calculate_gbox(lwg, &gbox_slow);
295  lwgeom_free(lwg);
296 #if 0
297  printf("\nCALC: %s\n", gbox_to_string(&gbox));
298  printf("GOOD: %s\n", gbox_to_string(&gbox_slow));
299  printf("line %d: diff %.9g\n", i, fabs(gbox.xmin - gbox_slow.xmin)+fabs(gbox.ymin - gbox_slow.ymin)+fabs(gbox.zmin - gbox_slow.zmin));
300  printf("------------\n");
301 #endif
302  CU_ASSERT_DOUBLE_EQUAL(gbox.xmin, gbox_slow.xmin, 0.00000001);
303  CU_ASSERT_DOUBLE_EQUAL(gbox.ymin, gbox_slow.ymin, 0.00000001);
304  CU_ASSERT_DOUBLE_EQUAL(gbox.zmin, gbox_slow.zmin, 0.00000001);
305  CU_ASSERT_DOUBLE_EQUAL(gbox.xmax, gbox_slow.xmax, 0.00000001);
306  CU_ASSERT_DOUBLE_EQUAL(gbox.ymax, gbox_slow.ymax, 0.00000001);
307  CU_ASSERT_DOUBLE_EQUAL(gbox.zmax, gbox_slow.zmax, 0.00000001);
308  }
309 
310 }
311 
312 /*
313 * Build LWGEOM on top of *aligned* structure so we can use the read-only
314 * point access methods on them.
315 static LWGEOM* lwgeom_over_gserialized(char *wkt)
316 {
317  LWGEOM *lwg;
318  GSERIALIZED *g;
319 
320  lwg = lwgeom_from_wkt(wkt, LW_PARSER_CHECK_NONE);
321  g = gserialized_from_lwgeom(lwg, 1, 0);
322  lwgeom_free(lwg);
323  return lwgeom_from_gserialized(g);
324 }
325 */
326 
327 static void edge_set(double lon1, double lat1, double lon2, double lat2, GEOGRAPHIC_EDGE *e)
328 {
329  e->start.lon = lon1;
330  e->start.lat = lat1;
331  e->end.lon = lon2;
332  e->end.lat = lat2;
333  edge_deg2rad(e);
334 }
335 
336 static void point_set(double lon, double lat, GEOGRAPHIC_POINT *p)
337 {
338  p->lon = lon;
339  p->lat = lat;
340  point_deg2rad(p);
341 }
342 
343 static void test_clairaut(void)
344 {
345 
346  GEOGRAPHIC_POINT gs, ge;
347  POINT3D vs, ve;
348  GEOGRAPHIC_POINT g_out_top, g_out_bottom, v_out_top, v_out_bottom;
349 
350  point_set(-45.0, 60.0, &gs);
351  point_set(135.0, 60.0, &ge);
352 
353  geog2cart(&gs, &vs);
354  geog2cart(&ge, &ve);
355 
356  clairaut_cartesian(&vs, &ve, &v_out_top, &v_out_bottom);
357  clairaut_geographic(&gs, &ge, &g_out_top, &g_out_bottom);
358 
359  CU_ASSERT_DOUBLE_EQUAL(v_out_top.lat, g_out_top.lat, 0.000001);
360  CU_ASSERT_DOUBLE_EQUAL(v_out_top.lon, g_out_top.lon, 0.000001);
361  CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lat, g_out_bottom.lat, 0.000001);
362  CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lon, g_out_bottom.lon, 0.000001);
363 
364  gs.lat = 1.3021240033804449;
365  ge.lat = 1.3021240033804449;
366  gs.lon = -1.3387392931438733;
367  ge.lon = 1.80285336044592;
368 
369  geog2cart(&gs, &vs);
370  geog2cart(&ge, &ve);
371 
372  clairaut_cartesian(&vs, &ve, &v_out_top, &v_out_bottom);
373  clairaut_geographic(&gs, &ge, &g_out_top, &g_out_bottom);
374 
375  CU_ASSERT_DOUBLE_EQUAL(v_out_top.lat, g_out_top.lat, 0.000001);
376  CU_ASSERT_DOUBLE_EQUAL(v_out_top.lon, g_out_top.lon, 0.000001);
377  CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lat, g_out_bottom.lat, 0.000001);
378  CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lon, g_out_bottom.lon, 0.000001);
379 }
380 
381 static void test_edge_intersection(void)
382 {
383  GEOGRAPHIC_EDGE e1, e2;
385  int rv;
386 
387  /* Covers case, end-to-end intersection */
388  edge_set(50, -10.999999999999998224, -10.0, 50.0, &e1);
389  edge_set(-10.0, 50.0, -10.272779983831613393, -16.937003313332997578, &e2);
390  rv = edge_intersection(&e1, &e2, &g);
391  CU_ASSERT_EQUAL(rv, LW_TRUE);
392 
393  /* Medford case, very short segment vs very long one */
394  e1.start.lat = 0.74123572595649878103;
395  e1.start.lon = -2.1496353191142714145;
396  e1.end.lat = 0.74123631950116664058;
397  e1.end.lon = -2.1496353248304860273;
398  e2.start.lat = 0.73856343764436815924;
399  e2.start.lon = -2.1461493501950630325;
400  e2.end.lat = 0.70971354024834598651;
401  e2.end.lon = 2.1082194552519770703;
402  rv = edge_intersection(&e1, &e2, &g);
403  CU_ASSERT_EQUAL(rv, LW_FALSE);
404 
405  /* Again, this time with a less exact input edge. */
406  edge_set(-123.165031277506, 42.4696787216231, -123.165031605021, 42.4697127292275, &e1);
407  rv = edge_intersection(&e1, &e2, &g);
408  CU_ASSERT_EQUAL(rv, LW_FALSE);
409 
410  /* Second Medford case, very short segment vs very long one
411  e1.start.lat = 0.73826546728290887156;
412  e1.start.lon = -2.14426380171833042;
413  e1.end.lat = 0.73826545883786642843;
414  e1.end.lon = -2.1442638997530165668;
415  e2.start.lat = 0.73775469118192538165;
416  e2.start.lon = -2.1436035534281718817;
417  e2.end.lat = 0.71021099548296817705;
418  e2.end.lon = 2.1065275171200439353;
419  rv = edge_intersection(e1, e2, &g);
420  CU_ASSERT_EQUAL(rv, LW_FALSE);
421  */
422 
423  /* Intersection at (0 0) */
424  edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
425  edge_set(0.0, -1.0, 0.0, 1.0, &e2);
426  rv = edge_intersection(&e1, &e2, &g);
427  point_rad2deg(&g);
428  CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
429  CU_ASSERT_DOUBLE_EQUAL(g.lon, 0.0, 0.00001);
430  CU_ASSERT_EQUAL(rv, LW_TRUE);
431 
432  /* No intersection at (0 0) */
433  edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
434  edge_set(0.0, -1.0, 0.0, -2.0, &e2);
435  rv = edge_intersection(&e1, &e2, &g);
436  CU_ASSERT_EQUAL(rv, LW_FALSE);
437 
438  /* End touches middle of segment at (0 0) */
439  edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
440  edge_set(0.0, -1.0, 0.0, 0.0, &e2);
441  rv = edge_intersection(&e1, &e2, &g);
442  point_rad2deg(&g);
443 #if 0
444  printf("\n");
445  printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon, e1.start.lat, e1.end.lon, e1.end.lat);
446  printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon, e2.start.lat, e2.end.lon, e2.end.lat);
447  printf("g = (%.15g %.15g)\n", g.lon, g.lat);
448  printf("rv = %d\n", rv);
449 #endif
450  CU_ASSERT_DOUBLE_EQUAL(g.lon, 0.0, 0.00001);
451  CU_ASSERT_EQUAL(rv, LW_TRUE);
452 
453  /* End touches end of segment at (0 0) */
454  edge_set(0.0, 0.0, 1.0, 0.0, &e1);
455  edge_set(0.0, -1.0, 0.0, 0.0, &e2);
456  rv = edge_intersection(&e1, &e2, &g);
457  point_rad2deg(&g);
458 #if 0
459  printf("\n");
460  printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon, e1.start.lat, e1.end.lon, e1.end.lat);
461  printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon, e2.start.lat, e2.end.lon, e2.end.lat);
462  printf("g = (%.15g %.15g)\n", g.lon, g.lat);
463  printf("rv = %d\n", rv);
464 #endif
465  CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
466  CU_ASSERT_DOUBLE_EQUAL(g.lon, 0.0, 0.00001);
467  CU_ASSERT_EQUAL(rv, LW_TRUE);
468 
469  /* Intersection at (180 0) */
470  edge_set(-179.0, -1.0, 179.0, 1.0, &e1);
471  edge_set(-179.0, 1.0, 179.0, -1.0, &e2);
472  rv = edge_intersection(&e1, &e2, &g);
473  point_rad2deg(&g);
474  CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
475  CU_ASSERT_DOUBLE_EQUAL(fabs(g.lon), 180.0, 0.00001);
476  CU_ASSERT_EQUAL(rv, LW_TRUE);
477 
478  /* Intersection at (180 0) */
479  edge_set(-170.0, 0.0, 170.0, 0.0, &e1);
480  edge_set(180.0, -10.0, 180.0, 10.0, &e2);
481  rv = edge_intersection(&e1, &e2, &g);
482  point_rad2deg(&g);
483  CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
484  CU_ASSERT_DOUBLE_EQUAL(fabs(g.lon), 180.0, 0.00001);
485  CU_ASSERT_EQUAL(rv, LW_TRUE);
486 
487  /* Intersection at north pole */
488  edge_set(-180.0, 80.0, 0.0, 80.0, &e1);
489  edge_set(90.0, 80.0, -90.0, 80.0, &e2);
490  rv = edge_intersection(&e1, &e2, &g);
491  point_rad2deg(&g);
492  CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
493  CU_ASSERT_EQUAL(rv, LW_TRUE);
494 
495  /* Equal edges return true */
496  edge_set(45.0, 10.0, 50.0, 20.0, &e1);
497  edge_set(45.0, 10.0, 50.0, 20.0, &e2);
498  rv = edge_intersection(&e1, &e2, &g);
499  point_rad2deg(&g);
500  CU_ASSERT_EQUAL(rv, LW_TRUE);
501 
502  /* Parallel edges (same great circle, different end points) return true */
503  edge_set(40.0, 0.0, 70.0, 0.0, &e1);
504  edge_set(60.0, 0.0, 50.0, 0.0, &e2);
505  rv = edge_intersection(&e1, &e2, &g);
506  point_rad2deg(&g);
507  CU_ASSERT_EQUAL(rv, 2); /* Hack, returning 2 as the 'co-linear' value */
508 
509  /* End touches arc at north pole */
510  edge_set(-180.0, 80.0, 0.0, 80.0, &e1);
511  edge_set(90.0, 80.0, -90.0, 90.0, &e2);
512  rv = edge_intersection(&e1, &e2, &g);
513  point_rad2deg(&g);
514 #if 0
515  printf("\n");
516  printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon, e1.start.lat, e1.end.lon, e1.end.lat);
517  printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon, e2.start.lat, e2.end.lon, e2.end.lat);
518  printf("g = (%.15g %.15g)\n", g.lon, g.lat);
519  printf("rv = %d\n", rv);
520 #endif
521  CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
522  CU_ASSERT_EQUAL(rv, LW_TRUE);
523 
524  /* End touches end at north pole */
525  edge_set(-180.0, 80.0, 0.0, 90.0, &e1);
526  edge_set(90.0, 80.0, -90.0, 90.0, &e2);
527  rv = edge_intersection(&e1, &e2, &g);
528  point_rad2deg(&g);
529 #if 0
530  printf("\n");
531  printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon, e1.start.lat, e1.end.lon, e1.end.lat);
532  printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon, e2.start.lat, e2.end.lon, e2.end.lat);
533  printf("g = (%.15g %.15g)\n", g.lon, g.lat);
534  printf("rv = %d\n", rv);
535 #endif
536  CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
537  CU_ASSERT_EQUAL(rv, LW_TRUE);
538 
539 }
540 
541 static void line2pts(const char *wkt, POINT3D *A1, POINT3D *A2)
542 {
544  POINTARRAY *pa;
545  POINT2D p1, p2;
546  GEOGRAPHIC_POINT g1, g2;
547  if ( ! l )
548  {
549  printf("BAD WKT FOUND in test_edge_intersects:\n %s\n\n", wkt);
550  exit(0);
551  }
552  pa = l->points;
553  getPoint2d_p(pa, 0, &p1);
554  getPoint2d_p(pa, 1, &p2);
555  geographic_point_init(p1.x, p1.y, &g1);
556  geographic_point_init(p2.x, p2.y, &g2);
557  geog2cart(&g1, A1);
558  geog2cart(&g2, A2);
559  lwline_free(l);
560  return;
561 }
562 
563 static void test_edge_intersects(void)
564 {
565  POINT3D A1, A2, B1, B2;
567  int rv;
568 
569  /* Covers case, end-to-end intersection */
570  line2pts("LINESTRING(50 -10.999999999999998224, -10.0 50.0)", &A1, &A2);
571  line2pts("LINESTRING(-10.0 50.0, -10.272779983831613393 -16.937003313332997578)", &B1, &B2);
572  rv = edge_intersects(&A1, &A2, &B1, &B2);
573  CU_ASSERT(rv & PIR_INTERSECTS);
574 
575  /* Medford case, very short segment vs very long one */
576  g.lat = 0.74123572595649878103;
577  g.lon = -2.1496353191142714145;
578  geog2cart(&g, &A1);
579  g.lat = 0.74123631950116664058;
580  g.lon = -2.1496353248304860273;
581  geog2cart(&g, &A2);
582  g.lat = 0.73856343764436815924;
583  g.lon = -2.1461493501950630325;
584  geog2cart(&g, &B1);
585  g.lat = 0.70971354024834598651;
586  g.lon = 2.1082194552519770703;
587  geog2cart(&g, &B2);
588  rv = edge_intersects(&A1, &A2, &B1, &B2);
589  CU_ASSERT(rv == 0);
590 
591  /* Second Medford case, very short segment vs very long one */
592  g.lat = 0.73826546728290887156;
593  g.lon = -2.14426380171833042;
594  geog2cart(&g, &A1);
595  g.lat = 0.73826545883786642843;
596  g.lon = -2.1442638997530165668;
597  geog2cart(&g, &A2);
598  g.lat = 0.73775469118192538165;
599  g.lon = -2.1436035534281718817;
600  geog2cart(&g, &B1);
601  g.lat = 0.71021099548296817705;
602  g.lon = 2.1065275171200439353;
603  geog2cart(&g, &B2);
604  rv = edge_intersects(&A1, &A2, &B1, &B2);
605  CU_ASSERT(rv == PIR_INTERSECTS);
606 
607  /* Again, this time with a less exact input edge. */
608  line2pts("LINESTRING(-123.165031277506 42.4696787216231, -123.165031605021 42.4697127292275)", &A1, &A2);
609  rv = edge_intersects(&A1, &A2, &B1, &B2);
610  CU_ASSERT(rv == 0);
611 
612  /* Intersection at (0 0) */
613  line2pts("LINESTRING(-1.0 0.0, 1.0 0.0)", &A1, &A2);
614  line2pts("LINESTRING(0.0 -1.0, 0.0 1.0)", &B1, &B2);
615  rv = edge_intersects(&A1, &A2, &B1, &B2);
616  CU_ASSERT(rv == PIR_INTERSECTS);
617 
618  /* No intersection at (0 0) */
619  line2pts("LINESTRING(-1.0 0.0, 1.0 0.0)", &A1, &A2);
620  line2pts("LINESTRING(0.0 -1.0, 0.0 -2.0)", &B1, &B2);
621  rv = edge_intersects(&A1, &A2, &B1, &B2);
622  CU_ASSERT(rv == 0);
623 
624  /* End touches middle of segment at (0 0) */
625  line2pts("LINESTRING(-1.0 0.0, 1.0 0.0)", &A1, &A2);
626  line2pts("LINESTRING(0.0 -1.0, 0.0 0.0)", &B1, &B2);
627  rv = edge_intersects(&A1, &A2, &B1, &B2);
628  CU_ASSERT(rv == (PIR_INTERSECTS|PIR_B_TOUCH_RIGHT) );
629 
630  /* End touches end of segment at (0 0) */
631  line2pts("LINESTRING(0.0 0.0, 1.0 0.0)", &A1, &A2);
632  line2pts("LINESTRING(0.0 -1.0, 0.0 0.0)", &B1, &B2);
633  rv = edge_intersects(&A1, &A2, &B1, &B2);
634  CU_ASSERT(rv == (PIR_INTERSECTS|PIR_B_TOUCH_RIGHT|PIR_A_TOUCH_RIGHT) );
635 
636  /* Intersection at (180 0) */
637  line2pts("LINESTRING(-179.0 -1.0, 179.0 1.0)", &A1, &A2);
638  line2pts("LINESTRING(-179.0 1.0, 179.0 -1.0)", &B1, &B2);
639  rv = edge_intersects(&A1, &A2, &B1, &B2);
640  CU_ASSERT(rv == PIR_INTERSECTS);
641 
642  /* Intersection at (180 0) */
643  line2pts("LINESTRING(-170.0 0.0, 170.0 0.0)", &A1, &A2);
644  line2pts("LINESTRING(180.0 -10.0, 180.0 10.0)", &B1, &B2);
645  rv = edge_intersects(&A1, &A2, &B1, &B2);
646  CU_ASSERT(rv == PIR_INTERSECTS);
647 
648  /* Intersection at north pole */
649  line2pts("LINESTRING(-180.0 80.0, 0.0 80.0)", &A1, &A2);
650  line2pts("LINESTRING(90.0 80.0, -90.0 80.0)", &B1, &B2);
651  rv = edge_intersects(&A1, &A2, &B1, &B2);
652  CU_ASSERT(rv == PIR_INTERSECTS);
653 
654  /* Equal edges return true */
655  line2pts("LINESTRING(45.0 10.0, 50.0 20.0)", &A1, &A2);
656  line2pts("LINESTRING(45.0 10.0, 50.0 20.0)", &B1, &B2);
657  rv = edge_intersects(&A1, &A2, &B1, &B2);
658  CU_ASSERT(rv & PIR_INTERSECTS);
659 
660  /* Parallel edges (same great circle, different end points) return true */
661  line2pts("LINESTRING(40.0 0.0, 70.0 0.0)", &A1, &A2);
662  line2pts("LINESTRING(60.0 0.0, 50.0 0.0)", &B1, &B2);
663  rv = edge_intersects(&A1, &A2, &B1, &B2);
664  CU_ASSERT(rv == (PIR_INTERSECTS|PIR_COLINEAR) );
665 
666  /* End touches arc at north pole */
667  line2pts("LINESTRING(-180.0 80.0, 0.0 80.0)", &A1, &A2);
668  line2pts("LINESTRING(90.0 80.0, -90.0 90.0)", &B1, &B2);
669  rv = edge_intersects(&A1, &A2, &B1, &B2);
670  CU_ASSERT(rv == (PIR_INTERSECTS|PIR_B_TOUCH_LEFT) );
671 
672  /* End touches end at north pole */
673  line2pts("LINESTRING(-180.0 80.0, 0.0 90.0)", &A1, &A2);
674  line2pts("LINESTRING(90.0 80.0, -90.0 90.0)", &B1, &B2);
675  rv = edge_intersects(&A1, &A2, &B1, &B2);
676  CU_ASSERT(rv == (PIR_INTERSECTS|PIR_B_TOUCH_LEFT|PIR_A_TOUCH_RIGHT) );
677 
678  /* Antipodal straddles. Great circles cross but at opposite */
679  /* sides of the globe */
680  /* #2534 */
681  /* http://www.gcmap.com/mapui?P=60N+90E-20S+90E%0D%0A0N+0E-90.04868865037885W+57.44011727050777S%0D%0A&MS=wls&DU=mi */
682  line2pts("LINESTRING(90.0 60.0, 90.0 -20.0)", &A1, &A2);
683  line2pts("LINESTRING(0.0 0.0, -90.04868865037885 -57.44011727050777)", &B1, &B2);
684  rv = edge_intersects(&A1, &A2, &B1, &B2);
685  CU_ASSERT(rv == 0);
686 
687  line2pts("LINESTRING(-5 0, 5 0)", &A1, &A2);
688  line2pts("LINESTRING(179 -5, 179 5)", &B1, &B2);
689  rv = edge_intersects(&A1, &A2, &B1, &B2);
690  CU_ASSERT(rv == 0);
691 
692  line2pts("LINESTRING(175 -85, 175 85)", &A1, &A2);
693  line2pts("LINESTRING(65 0, -105 0)", &B1, &B2);
694  rv = edge_intersects(&A1, &A2, &B1, &B2);
695  CU_ASSERT(rv == 0);
696 
697  line2pts("LINESTRING(175 -85, 175 85)", &A1, &A2);
698  line2pts("LINESTRING(45 0, -125 0)", &B1, &B2);
699  rv = edge_intersects(&A1, &A2, &B1, &B2);
700  CU_ASSERT(rv == 0);
701 
702 }
703 
705 {
706  GEOGRAPHIC_EDGE e;
708  GEOGRAPHIC_POINT closest;
709  double d;
710 
711  /* closest point at origin, one degree away */
712  edge_set(-50.0, 0.0, 50.0, 0.0, &e);
713  point_set(0.0, 1.0, &g);
714  d = edge_distance_to_point(&e, &g, 0);
715  CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001);
716 
717  /* closest point at origin, one degree away */
718  edge_set(-50.0, 0.0, 50.0, 0.0, &e);
719  point_set(0.0, 2.0, &g);
720  d = edge_distance_to_point(&e, &g, &closest);
721 #if 0
722  printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e.start.lon, e.start.lat, e.end.lon, e.end.lat);
723  printf("POINT(%.9g %.9g)\n", g.lon, g.lat);
724  printf("\nDISTANCE == %.8g\n", d);
725 #endif
726  CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 90.0, 0.00001);
727  CU_ASSERT_DOUBLE_EQUAL(closest.lat, 0.0, 0.00001);
728  CU_ASSERT_DOUBLE_EQUAL(closest.lon, 0.0, 0.00001);
729 
730  /* Ticket #2351 */
731  edge_set(149.386990599235, -26.3567415843982, 149.386990599247, -26.3567415843965, &e);
732  point_set(149.386990599235, -26.3567415843982, &g);
733  d = edge_distance_to_point(&e, &g, &closest);
734  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
735  // printf("CLOSE POINT(%g %g)\n", closest.lon, closest.lat);
736  // printf(" ORIG POINT(%g %g)\n", g.lon, g.lat);
737  CU_ASSERT_DOUBLE_EQUAL(g.lat, closest.lat, 0.00001);
738  CU_ASSERT_DOUBLE_EQUAL(g.lon, closest.lon, 0.00001);
739 }
740 
741 static void test_edge_distance_to_edge(void)
742 {
743  GEOGRAPHIC_EDGE e1, e2;
744  GEOGRAPHIC_POINT c1, c2;
745  double d;
746 
747  /* closest point at origin, one degree away */
748  edge_set(-50.0, 0.0, 50.0, 0.0, &e1);
749  edge_set(-5.0, 20.0, 0.0, 1.0, &e2);
750  d = edge_distance_to_edge(&e1, &e2, &c1, &c2);
751 #if 0
752  printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e1.start.lon, e1.start.lat, e1.end.lon, e1.end.lat);
753  printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e2.start.lon, e2.start.lat, e2.end.lon, e2.end.lat);
754  printf("\nDISTANCE == %.8g\n", d);
755 #endif
756  CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001);
757  CU_ASSERT_DOUBLE_EQUAL(c1.lat, 0.0, 0.00001);
758  CU_ASSERT_DOUBLE_EQUAL(c2.lat, M_PI / 180.0, 0.00001);
759  CU_ASSERT_DOUBLE_EQUAL(c1.lon, 0.0, 0.00001);
760  CU_ASSERT_DOUBLE_EQUAL(c2.lon, 0.0, 0.00001);
761 }
762 
763 
764 
765 /*
766 * Build LWGEOM on top of *aligned* structure so we can use the read-only
767 * point access methods on them.
768 */
770 {
771  LWGEOM *lwg;
772 
774  FLAGS_SET_GEODETIC(lwg->flags, 1);
775  *g = gserialized_from_lwgeom(lwg, 1, 0);
776  lwgeom_free(lwg);
777  return lwgeom_from_gserialized(*g);
778 }
779 
780 static void test_lwgeom_check_geodetic(void)
781 {
782  LWGEOM *geom;
783  int i = 0;
784 
785  char ewkt[][512] =
786  {
787  "POINT(0 0.2)",
788  "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
789  "SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
790  "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
791  "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
792  "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
793  "POINT(0 220.2)",
794  "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
795  "SRID=1;MULTILINESTRING((-1 -131,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
796  "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
797  "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,211 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
798  "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1111 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
799  };
800 
801  for ( i = 0; i < 6; i++ )
802  {
803  GSERIALIZED *g;
804  geom = lwgeom_over_gserialized(ewkt[i], &g);
805  CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), LW_TRUE);
806  lwgeom_free(geom);
807  lwfree(g);
808  }
809 
810  for ( i = 6; i < 12; i++ )
811  {
812  GSERIALIZED *g;
813  //char *out_ewkt;
814  geom = lwgeom_over_gserialized(ewkt[i], &g);
815  CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), LW_FALSE);
816  //out_ewkt = lwgeom_to_ewkt(geom);
817  //printf("%s\n", out_ewkt);
818  lwgeom_free(geom);
819  lwfree(g);
820  }
821 
822 }
823 
824 /*
825 static void test_gbox_calculation(void)
826 {
827 
828  LWGEOM *geom;
829  int i = 0;
830  GBOX *gbox = gbox_new(gflags(0,0,0));
831  BOX3D *box3d;
832 
833  char ewkt[][512] =
834  {
835  "POINT(0 0.2)",
836  "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
837  "SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
838  "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
839  "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
840  "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
841  "POINT(0 220.2)",
842  "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
843  "SRID=1;MULTILINESTRING((-1 -131,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
844  "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
845  "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,211 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
846  "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1111 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
847  };
848 
849  for ( i = 0; i < 6; i++ )
850  {
851  GSERIALIZED *g;
852  geom = lwgeom_over_gserialized(ewkt[i], &g);
853  lwgeom_calculate_gbox_cartesian(geom, gbox);
854  box3d = lwgeom_compute_box3d(geom);
855  //printf("%g %g\n", gbox->xmin, box3d->xmin);
856  CU_ASSERT_EQUAL(gbox->xmin, box3d->xmin);
857  CU_ASSERT_EQUAL(gbox->xmax, box3d->xmax);
858  CU_ASSERT_EQUAL(gbox->ymin, box3d->ymin);
859  CU_ASSERT_EQUAL(gbox->ymax, box3d->ymax);
860  lwgeom_free(geom);
861  lwfree(box3d);
862  lwfree(g);
863  }
864  lwfree(gbox);
865 }
866 */
867 
869 {
870  LWGEOM *geom;
871  GSERIALIZED *g;
872  uint32_t type;
873  double *inspect; /* To poke right into the blob. */
874 
875  geom = lwgeom_from_wkt("POINT(0 0.2)", LW_PARSER_CHECK_NONE);
876  FLAGS_SET_GEODETIC(geom->flags, 1);
877  g = gserialized_from_lwgeom(geom, 1, 0);
878  type = gserialized_get_type(g);
879  CU_ASSERT_EQUAL( type, POINTTYPE );
880  inspect = (double*)g;
881  CU_ASSERT_EQUAL(inspect[3], 0.2);
882  lwgeom_free(geom);
883  lwfree(g);
884 
885  geom = lwgeom_from_wkt("POLYGON((-1 -1, -1 2.5, 2 2, 2 -1, -1 -1), (0 0, 0 1, 1 1, 1 0, 0 0))", LW_PARSER_CHECK_NONE);
886  FLAGS_SET_GEODETIC(geom->flags, 1);
887  g = gserialized_from_lwgeom(geom, 1, 0);
888  type = gserialized_get_type(g);
889  CU_ASSERT_EQUAL( type, POLYGONTYPE );
890  inspect = (double*)g;
891  CU_ASSERT_EQUAL(inspect[9], 2.5);
892  lwgeom_free(geom);
893  lwfree(g);
894 
895  geom = lwgeom_from_wkt("MULTILINESTRING((0 0, 1 1),(0 0.1, 1 1))", LW_PARSER_CHECK_NONE);
896  FLAGS_SET_GEODETIC(geom->flags, 1);
897  g = gserialized_from_lwgeom(geom, 1, 0);
898  type = gserialized_get_type(g);
899  CU_ASSERT_EQUAL( type, MULTILINETYPE );
900  inspect = (double*)g;
901  CU_ASSERT_EQUAL(inspect[12], 0.1);
902  lwgeom_free(geom);
903  lwfree(g);
904 
905 }
906 
908 {
909  LWGEOM *lwg;
910  LWPOLY *poly;
911  POINT2D pt_to_test;
912  POINT2D pt_outside;
913  int result;
914 
915  /* Small polygon and huge distance between outside point and close-but-not-quite-inside point. Should return LW_FALSE. Pretty degenerate case. */
916  lwg = lwgeom_from_hexwkb("0103000020E61000000100000025000000ACAD6F91DDB65EC03F84A86D57264540CCABC279DDB65EC0FCE6926B57264540B6DEAA62DDB65EC0A79F6B63572645402E0BE84CDDB65EC065677155572645405D0B1D39DDB65EC0316310425726454082B5DB27DDB65EC060A4E12957264540798BB619DDB65EC0C393A10D57264540D4BC160FDDB65EC0BD0320EE56264540D7AC4E08DDB65EC096C862CC56264540AFD29205DDB65EC02A1F68A956264540363AFA06DDB65EC0722E418656264540B63A780CDDB65EC06E9B0064562645409614E215DDB65EC0E09DA84356264540FF71EF22DDB65EC0B48145265626454036033F33DDB65EC081B8A60C5626454066FB4546DDB65EC08A47A6F7552645409061785BDDB65EC0F05AE0E755264540D4B63772DDB65EC05C86CEDD55264540D2E4C689DDB65EC09B6EBFD95526454082E573A1DDB65EC0C90BD5DB552645401ABE85B8DDB65EC06692FCE35526454039844ECEDDB65EC04D8AF6F155264540928319E2DDB65EC0AD8D570556264540D31055F3DDB65EC02D618F1D56264540343B7A01DEB65EC0EB70CF3956264540920A1A0CDEB65EC03B00515956264540911BE212DEB65EC0E43A0E7B56264540E3F69D15DEB65EC017E4089E562645408D903614DEB65EC0F0D42FC1562645402191B80EDEB65EC0586870E35626454012B84E05DEB65EC09166C80357264540215B41F8DDB65EC08F832B21572645408392F7E7DDB65EC01138C13A57264540F999F0D4DDB65EC0E4A9C14F57264540AC3FB8BFDDB65EC0EED6875F57264540D3DCFEA8DDB65EC04F6C996957264540ACAD6F91DDB65EC03F84A86D57264540", LW_PARSER_CHECK_NONE);
917  poly = (LWPOLY*)lwg;
918  pt_to_test.x = -122.819436560680316;
919  pt_to_test.y = 42.2702301207017328;
920  pt_outside.x = 120.695136159150778;
921  pt_outside.y = 40.6920926049588516;
922  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
923  CU_ASSERT_EQUAL(result, LW_FALSE);
924  lwgeom_free(lwg);
925 
926  /* Point on ring between vertexes case */
927  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
928  poly = (LWPOLY*)lwg;
929  pt_to_test.x = 1.1;
930  pt_to_test.y = 1.05;
931  pt_outside.x = 1.2;
932  pt_outside.y = 1.05;
933  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
934  CU_ASSERT_EQUAL(result, LW_TRUE);
935  lwgeom_free(lwg);
936 
937  /* Simple containment case */
938  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
939  poly = (LWPOLY*)lwg;
940  pt_to_test.x = 1.05;
941  pt_to_test.y = 1.05;
942  pt_outside.x = 1.2;
943  pt_outside.y = 1.15;
944  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
945  CU_ASSERT_EQUAL(result, LW_TRUE);
946  lwgeom_free(lwg);
947 
948  /* Less Simple containment case. */
949  /* Interior point quite close to boundary and stab line going through bottom edge vertex */
950  /* This breaks the "extend-it" trick of handling vertex crossings */
951  /* It should also break the "lowest end" trick. */
952  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.05 0.95, 1.0 1.0))", LW_PARSER_CHECK_NONE);
953  poly = (LWPOLY*)lwg;
954  pt_to_test.x = 1.05;
955  pt_to_test.y = 1.00;
956  pt_outside.x = 1.05;
957  pt_outside.y = 0.5;
958  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
959  CU_ASSERT_EQUAL(result, LW_TRUE);
960  lwgeom_free(lwg);
961 
962  /* Simple noncontainment case */
963  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
964  poly = (LWPOLY*)lwg;
965  pt_to_test.x = 1.05;
966  pt_to_test.y = 1.15;
967  pt_outside.x = 1.2;
968  pt_outside.y = 1.2;
969  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
970  CU_ASSERT_EQUAL(result, LW_FALSE);
971  lwgeom_free(lwg);
972 
973  /* Harder noncontainment case */
974  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
975  poly = (LWPOLY*)lwg;
976  pt_to_test.x = 1.05;
977  pt_to_test.y = 0.9;
978  pt_outside.x = 1.2;
979  pt_outside.y = 1.05;
980  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
981  CU_ASSERT_EQUAL(result, LW_FALSE);
982  lwgeom_free(lwg);
983 
984  /* Harder containment case */
985  lwg = lwgeom_from_wkt("POLYGON((0 0, 0 2, 1 2, 0 3, 2 3, 0 4, 3 5, 0 6, 6 10, 6 1, 0 0))", LW_PARSER_CHECK_NONE);
986  poly = (LWPOLY*)lwg;
987  pt_to_test.x = 1.0;
988  pt_to_test.y = 1.0;
989  pt_outside.x = 1.0;
990  pt_outside.y = 10.0;
991  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
992  CU_ASSERT_EQUAL(result, LW_TRUE);
993  lwgeom_free(lwg);
994 
995  /* Point on ring at first vertex case */
996  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
997  poly = (LWPOLY*)lwg;
998  pt_to_test.x = 1.0;
999  pt_to_test.y = 1.0;
1000  pt_outside.x = 1.2;
1001  pt_outside.y = 1.05;
1002  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1003  CU_ASSERT_EQUAL(result, LW_TRUE);
1004  lwgeom_free(lwg);
1005 
1006  /* Point on ring at vertex case */
1007  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
1008  poly = (LWPOLY*)lwg;
1009  pt_to_test.x = 1.0;
1010  pt_to_test.y = 1.1;
1011  pt_outside.x = 1.2;
1012  pt_outside.y = 1.05;
1013  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1014  CU_ASSERT_EQUAL(result, LW_TRUE);
1015  lwgeom_free(lwg);
1016 
1017  /* Co-linear crossing case for point-in-polygon test, should return LW_TRUE */
1018  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.2, 1.2 1.2, 1.2 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
1019  poly = (LWPOLY*)lwg;
1020  pt_to_test.x = 1.1;
1021  pt_to_test.y = 1.05;
1022  pt_outside.x = 1.1;
1023  pt_outside.y = 1.3;
1024  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1025  CU_ASSERT_EQUAL(result, LW_TRUE);
1026  lwgeom_free(lwg);
1027 
1028  /* Co-linear grazing case for point-in-polygon test, should return LW_FALSE */
1029  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.2, 1.2 1.2, 1.2 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
1030  poly = (LWPOLY*)lwg;
1031  pt_to_test.x = 1.0;
1032  pt_to_test.y = 0.0;
1033  pt_outside.x = 1.0;
1034  pt_outside.y = 2.0;
1035  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1036  CU_ASSERT_EQUAL(result, LW_FALSE);
1037  lwgeom_free(lwg);
1038 
1039  /* Grazing case for point-in-polygon test, should return LW_FALSE */
1040  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 2.0, 1.5 1.5, 1.0 1.0))", LW_PARSER_CHECK_NONE);
1041  poly = (LWPOLY*)lwg;
1042  pt_to_test.x = 1.5;
1043  pt_to_test.y = 1.0;
1044  pt_outside.x = 1.5;
1045  pt_outside.y = 2.0;
1046  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1047  CU_ASSERT_EQUAL(result, LW_FALSE);
1048  lwgeom_free(lwg);
1049 
1050  /* Grazing case at first point for point-in-polygon test, should return LW_FALSE */
1051  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 2.0 3.0, 2.0 0.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
1052  poly = (LWPOLY*)lwg;
1053  pt_to_test.x = 1.0;
1054  pt_to_test.y = 0.0;
1055  pt_outside.x = 1.0;
1056  pt_outside.y = 2.0;
1057  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1058  CU_ASSERT_EQUAL(result, LW_FALSE);
1059  lwgeom_free(lwg);
1060 
1061  /* Outside multi-crossing case for point-in-polygon test, should return LW_FALSE */
1062  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.2, 1.2 1.2, 1.2 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
1063  poly = (LWPOLY*)lwg;
1064  pt_to_test.x = 0.99;
1065  pt_to_test.y = 0.99;
1066  pt_outside.x = 1.21;
1067  pt_outside.y = 1.21;
1068  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1069  CU_ASSERT_EQUAL(result, LW_FALSE);
1070  lwgeom_free(lwg);
1071 
1072  /* Inside multi-crossing case for point-in-polygon test, should return LW_TRUE */
1073  lwg = lwgeom_from_wkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.2, 1.2 1.2, 1.2 1.0, 1.0 1.0))", LW_PARSER_CHECK_NONE);
1074  poly = (LWPOLY*)lwg;
1075  pt_to_test.x = 1.11;
1076  pt_to_test.y = 1.11;
1077  pt_outside.x = 1.21;
1078  pt_outside.y = 1.21;
1079  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1080  CU_ASSERT_EQUAL(result, LW_TRUE);
1081  lwgeom_free(lwg);
1082 
1083  /* Point on vertex of ring */
1084  lwg = lwgeom_from_wkt("POLYGON((-9 50,51 -11,-10 50,-9 50))", LW_PARSER_CHECK_NONE);
1085  poly = (LWPOLY*)lwg;
1086  pt_to_test.x = -10.0;
1087  pt_to_test.y = 50.0;
1088  pt_outside.x = -10.2727799838316134;
1089  pt_outside.y = -16.9370033133329976;
1090  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1091  CU_ASSERT_EQUAL(result, LW_TRUE);
1092  lwgeom_free(lwg);
1093 
1094 }
1095 
1097 {
1098  LWPOLY *poly;
1099  LWGEOM *lwg;
1100  POINT2D pt_to_test;
1101  int result;
1102 
1103  lwg = lwgeom_from_wkt("POLYGON((-9 50,51 -11,-10 50,-9 50))", LW_PARSER_CHECK_NONE);
1104  poly = (LWPOLY*)lwg;
1105  pt_to_test.x = -10.0;
1106  pt_to_test.y = 50.0;
1107  result = lwpoly_covers_point2d(poly, &pt_to_test);
1108  CU_ASSERT_EQUAL(result, LW_TRUE);
1109  lwgeom_free(lwg);
1110 
1111  /* Great big ring */
1112  lwg = lwgeom_from_wkt("POLYGON((-40.0 52.0, 102.0 -6.0, -67.0 -29.0, -40.0 52.0))", LW_PARSER_CHECK_NONE);
1113  poly = (LWPOLY*)lwg;
1114  pt_to_test.x = 4.0;
1115  pt_to_test.y = 11.0;
1116  result = lwpoly_covers_point2d(poly, &pt_to_test);
1117  CU_ASSERT_EQUAL(result, LW_TRUE);
1118  lwgeom_free(lwg);
1119 
1120 }
1121 
1123 {
1125  LWPOLY *poly = (LWPOLY*)lwg;
1126  POINTARRAY *pa = poly->rings[0];
1127  POINT2D pt_outside, pt_to_test;
1128  int rv;
1129 
1130  pt_to_test.x = -95.900000000000006;
1131  pt_to_test.y = 42.899999999999999;
1132  pt_outside.x = -96.381873780830645;
1133  pt_outside.y = 40.185394449416371;
1134 
1135  rv = ptarray_contains_point_sphere(pa, &pt_outside, &pt_to_test);
1136  CU_ASSERT_EQUAL(rv, LW_TRUE);
1137 
1138  lwgeom_free(lwg);
1139 }
1140 
1141 
1143 {
1144  LWGEOM *lwg1, *lwg2;
1145  double d;
1146  SPHEROID s;
1147 
1148  /* Init and force spherical */
1149  spheroid_init(&s, 6378137.0, 6356752.314245179498);
1150  s.a = s.b = s.radius;
1151 
1152  /* Line/line distance, 1 degree apart */
1153  lwg1 = lwgeom_from_wkt("LINESTRING(-30 10, -20 5, -10 3, 0 1)", LW_PARSER_CHECK_NONE);
1154  lwg2 = lwgeom_from_wkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", LW_PARSER_CHECK_NONE);
1155  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1156  CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001);
1157  lwgeom_free(lwg1);
1158  lwgeom_free(lwg2);
1159 
1160  /* Line/line distance, crossing, 0.0 apart */
1161  lwg1 = lwgeom_from_wkt("LINESTRING(-30 10, -20 5, -10 3, 0 1)", LW_PARSER_CHECK_NONE);
1162  lwg2 = lwgeom_from_wkt("LINESTRING(-10 -5, -5 20, 5 0, 10 -5)", LW_PARSER_CHECK_NONE);
1163  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1164  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1165  lwgeom_free(lwg1);
1166  lwgeom_free(lwg2);
1167 
1168  /* Line/point distance, 1 degree apart */
1169  lwg1 = lwgeom_from_wkt("POINT(-4 1)", LW_PARSER_CHECK_NONE);
1170  lwg2 = lwgeom_from_wkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", LW_PARSER_CHECK_NONE);
1171  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1172  CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001);
1173  lwgeom_free(lwg1);
1174  lwgeom_free(lwg2);
1175 
1176  lwg1 = lwgeom_from_wkt("POINT(-4 1)", LW_PARSER_CHECK_NONE);
1177  lwg2 = lwgeom_from_wkt("POINT(-4 -1)", LW_PARSER_CHECK_NONE);
1178  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1179  CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 90.0, 0.00001);
1180  lwgeom_free(lwg1);
1181  lwgeom_free(lwg2);
1182 
1183  /* Poly/point distance, point inside polygon, 0.0 apart */
1184  lwg1 = lwgeom_from_wkt("POLYGON((-4 1, -3 5, 1 2, 1.5 -5, -4 1))", LW_PARSER_CHECK_NONE);
1185  lwg2 = lwgeom_from_wkt("POINT(-1 -1)", LW_PARSER_CHECK_NONE);
1186  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1187  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1188  lwgeom_free(lwg1);
1189  lwgeom_free(lwg2);
1190 
1191  /* Poly/point distance, point inside polygon hole, 1 degree apart */
1192  lwg1 = lwgeom_from_wkt("POLYGON((-4 -4, -4 4, 4 4, 4 -4, -4 -4), (-2 -2, -2 2, 2 2, 2 -2, -2 -2))", LW_PARSER_CHECK_NONE);
1193  lwg2 = lwgeom_from_wkt("POINT(-1 -1)", LW_PARSER_CHECK_NONE);
1194  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1195  CU_ASSERT_DOUBLE_EQUAL(d, 111178.142466, 0.1);
1196  lwgeom_free(lwg1);
1197  lwgeom_free(lwg2);
1198 
1199  /* Poly/point distance, point on hole boundary, 0.0 apart */
1200  lwg1 = lwgeom_from_wkt("POLYGON((-4 -4, -4 4, 4 4, 4 -4, -4 -4), (-2 -2, -2 2, 2 2, 2 -2, -2 -2))", LW_PARSER_CHECK_NONE);
1201  lwg2 = lwgeom_from_wkt("POINT(2 2)", LW_PARSER_CHECK_NONE);
1202  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1203  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1204  lwgeom_free(lwg1);
1205  lwgeom_free(lwg2);
1206 
1207  /* Medford test case #1 */
1208  lwg1 = lwgeom_from_hexwkb("0105000020E610000001000000010200000002000000EF7B8779C7BD5EC0FD20D94B852845400E539C62B9BD5EC0F0A5BE767C284540", LW_PARSER_CHECK_NONE);
1209  lwg2 = lwgeom_from_hexwkb("0106000020E61000000100000001030000000100000007000000280EC3FB8CCA5EC0A5CDC747233C45402787C8F58CCA5EC0659EA2761E3C45400CED58DF8FCA5EC0C37FAE6E1E3C4540AE97B8E08FCA5EC00346F58B1F3C4540250359FD8ECA5EC05460628E1F3C45403738F4018FCA5EC05DC84042233C4540280EC3FB8CCA5EC0A5CDC747233C4540", LW_PARSER_CHECK_NONE);
1210  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1211  CU_ASSERT_DOUBLE_EQUAL(d, 23630.8003, 0.1);
1212  lwgeom_free(lwg1);
1213  lwgeom_free(lwg2);
1214 
1215  /* Ticket #2351 */
1216  lwg1 = lwgeom_from_wkt("LINESTRING(149.386990599235 -26.3567415843982,149.386990599247 -26.3567415843965)", LW_PARSER_CHECK_NONE);
1217  lwg2 = lwgeom_from_wkt("POINT(149.386990599235 -26.3567415843982)", LW_PARSER_CHECK_NONE);
1218  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1219  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1220  lwgeom_free(lwg1);
1221  lwgeom_free(lwg2);
1222 
1223  /* Ticket #2638, no "M" */
1224  lwg1 = lwgeom_from_wkt("LINESTRING (-41.0821 50.3036,50 -41)", LW_PARSER_CHECK_NONE);
1225  lwg2 = lwgeom_from_wkt("POLYGON((0 0,10 0,10 10,0 10,0 0),(5 5,7 5,7 7,5 7,5 5))", LW_PARSER_CHECK_NONE);
1226  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1227  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1228  lwgeom_free(lwg1);
1229  lwgeom_free(lwg2);
1230 
1231  /* Ticket #2638, with "M" */
1232  lwg1 = lwgeom_from_wkt("LINESTRING M (-41.0821 50.3036 1,50 -41 1)", LW_PARSER_CHECK_NONE);
1233  lwg2 = lwgeom_from_wkt("POLYGON M ((0 0 2,10 0 1,10 10 -2,0 10 -5,0 0 -5),(5 5 6,7 5 6,7 7 6,5 7 10,5 5 -2))", LW_PARSER_CHECK_NONE);
1234  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1235  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1236  lwgeom_free(lwg1);
1237  lwgeom_free(lwg2);
1238 }
1239 
1240 static void test_spheroid_distance(void)
1241 {
1242  GEOGRAPHIC_POINT g1, g2;
1243  double d;
1244  SPHEROID s;
1245 
1246  /* Init to WGS84 */
1247  spheroid_init(&s, 6378137.0, 6356752.314245179498);
1248 
1249  /* One vertical degree */
1250  point_set(0.0, 0.0, &g1);
1251  point_set(0.0, 1.0, &g2);
1252  d = spheroid_distance(&g1, &g2, &s);
1253  CU_ASSERT_DOUBLE_EQUAL(d, 110574.388615329, 0.001);
1254 
1255  /* Ten horizontal degrees */
1256  point_set(-10.0, 0.0, &g1);
1257  point_set(0.0, 0.0, &g2);
1258  d = spheroid_distance(&g1, &g2, &s);
1259  CU_ASSERT_DOUBLE_EQUAL(d, 1113194.90793274, 0.001);
1260 
1261  /* One horizonal degree */
1262  point_set(-1.0, 0.0, &g1);
1263  point_set(0.0, 0.0, &g2);
1264  d = spheroid_distance(&g1, &g2, &s);
1265  CU_ASSERT_DOUBLE_EQUAL(d, 111319.490779, 0.001);
1266 
1267  /* Around world w/ slight bend */
1268  point_set(-180.0, 0.0, &g1);
1269  point_set(0.0, 1.0, &g2);
1270  d = spheroid_distance(&g1, &g2, &s);
1271  CU_ASSERT_DOUBLE_EQUAL(d, 19893357.0704483, 0.001);
1272 
1273  /* Up to pole */
1274  point_set(-180.0, 0.0, &g1);
1275  point_set(0.0, 90.0, &g2);
1276  d = spheroid_distance(&g1, &g2, &s);
1277  CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7295318, 0.001);
1278 
1279 }
1280 
1281 static void test_spheroid_area(void)
1282 {
1283  LWGEOM *lwg;
1284  GBOX gbox;
1285  double a1, a2;
1286  SPHEROID s;
1287 
1288  /* Init to WGS84 */
1290 
1291  gbox.flags = gflags(0, 0, 1);
1292 
1293  /* Medford lot test polygon */
1294  lwg = lwgeom_from_wkt("POLYGON((-122.848227067007 42.5007249610493,-122.848309475585 42.5007179884263,-122.848327688675 42.500835880696,-122.848245279942 42.5008428533324,-122.848227067007 42.5007249610493))", LW_PARSER_CHECK_NONE);
1295  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1296  a1 = lwgeom_area_sphere(lwg, &s);
1297  a2 = lwgeom_area_spheroid(lwg, &s);
1298  //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
1299  CU_ASSERT_DOUBLE_EQUAL(a1, 89.7127703297, 0.1); /* sphere */
1300  CU_ASSERT_DOUBLE_EQUAL(a2, 89.8684316032, 0.1); /* spheroid */
1301  lwgeom_free(lwg);
1302 
1303  /* Big-ass polygon */
1304  lwg = lwgeom_from_wkt("POLYGON((-2 3, -2 4, -1 4, -1 3, -2 3))", LW_PARSER_CHECK_NONE);
1305  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1306  a1 = lwgeom_area_sphere(lwg, &s);
1307  a2 = lwgeom_area_spheroid(lwg, &s);
1308  //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
1309  CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */
1310  CU_ASSERT_DOUBLE_EQUAL(a2, 12286574431.9, 10.0); /* spheroid */
1311  lwgeom_free(lwg);
1312 
1313  /* One-degree square */
1314  lwg = lwgeom_from_wkt("POLYGON((8.5 2,8.5 1,9.5 1,9.5 2,8.5 2))", LW_PARSER_CHECK_NONE);
1315  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1316  a1 = lwgeom_area_sphere(lwg, &s);
1317  a2 = lwgeom_area_spheroid(lwg, &s);
1318  //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
1319  CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */
1320  CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
1321  lwgeom_free(lwg);
1322 
1323  /* One-degree square *near* dateline */
1324  lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,178.5 1,178.5 2,179.5 2))", LW_PARSER_CHECK_NONE);
1325  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1326  a1 = lwgeom_area_sphere(lwg, &s);
1327  a2 = lwgeom_area_spheroid(lwg, &s);
1328  //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
1329  CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */
1330  CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
1331  lwgeom_free(lwg);
1332 
1333  /* One-degree square *across* dateline */
1334  lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", LW_PARSER_CHECK_NONE);
1335  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1336  a1 = lwgeom_area_sphere(lwg, &s);
1337  a2 = lwgeom_area_spheroid(lwg, &s);
1338  //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
1339  CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.3679, 10.0); /* sphere */
1340  CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
1341  lwgeom_free(lwg);
1342 }
1343 
1344 static void test_gbox_utils(void)
1345 {
1346  LWGEOM *lwg;
1347  GBOX gbox;
1348  double a1, a2;
1349  SPHEROID s;
1350  POINT2D pt;
1351 
1352  /* Init to WGS84 */
1354 
1355  gbox.flags = gflags(0, 0, 1);
1356 
1357  /* One-degree square by equator */
1358  lwg = lwgeom_from_wkt("POLYGON((1 20,1 21,2 21,2 20,1 20))", LW_PARSER_CHECK_NONE);
1359  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1360  a1 = gbox_angular_width(&gbox);
1361  a2 = gbox_angular_height(&gbox);
1362  CU_ASSERT_DOUBLE_EQUAL(a1, 0.0177951, 0.0000001);
1363  CU_ASSERT_DOUBLE_EQUAL(a2, 0.017764, 0.0000001);
1364  lwgeom_free(lwg);
1365 
1366  /* One-degree square *across* dateline */
1367  lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", LW_PARSER_CHECK_NONE);
1368  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1369  a1 = gbox_angular_width(&gbox);
1370  a2 = gbox_angular_height(&gbox);
1371  //printf("a1=%g a2=%g\n", a1, a2);
1372  CU_ASSERT_DOUBLE_EQUAL(a1, 0.0174613, 0.0000001);
1373  CU_ASSERT_DOUBLE_EQUAL(a2, 0.0174553, 0.0000001);
1374  lwgeom_free(lwg);
1375 
1376  /* One-degree square *across* dateline */
1377  lwg = lwgeom_from_wkt("POLYGON((178.5 2,178.5 1,-179.5 1,-179.5 2,178.5 2))", LW_PARSER_CHECK_NONE);
1378  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1379  gbox_centroid(&gbox, &pt);
1380  //printf("POINT(%g %g)\n", pt.x, pt.y);
1381  CU_ASSERT_DOUBLE_EQUAL(pt.x, 179.5, 0.0001);
1382  CU_ASSERT_DOUBLE_EQUAL(pt.y, 1.50024, 0.0001);
1383  lwgeom_free(lwg);
1384 
1385 }
1386 
1387 static void test_vector_angle(void)
1388 {
1389  POINT3D p1, p2;
1390  double angle;
1391 
1392  memset(&p1, 0, sizeof(POINT3D));
1393  memset(&p2, 0, sizeof(POINT3D));
1394 
1395  p1.x = 1.0;
1396  p2.y = 1.0;
1397  angle = vector_angle(&p1, &p2);
1398  CU_ASSERT_DOUBLE_EQUAL(angle, M_PI/2, 0.00001);
1399 
1400  p1.x = p2.y = 0.0;
1401  p1.y = 1.0;
1402  p2.x = 1.0;
1403  angle = vector_angle(&p1, &p2);
1404  CU_ASSERT_DOUBLE_EQUAL(angle, M_PI/2, 0.00001);
1405 
1406  p2.y = p2.x = 1.0;
1407  normalize(&p2);
1408  angle = vector_angle(&p1, &p2);
1409  CU_ASSERT_DOUBLE_EQUAL(angle, M_PI/4, 0.00001);
1410 
1411  p2.x = p2.y = p2.z = 1.0;
1412  normalize(&p2);
1413  angle = vector_angle(&p1, &p2);
1414  CU_ASSERT_DOUBLE_EQUAL(angle, 0.955317, 0.00001);
1415  //printf ("angle = %g\n\n", angle);
1416 }
1417 
1418 static void test_vector_rotate(void)
1419 {
1420  POINT3D p1, p2, n;
1421  double angle;
1422 
1423  memset(&p1, 0, sizeof(POINT3D));
1424  memset(&p2, 0, sizeof(POINT3D));
1425  memset(&n, 0, sizeof(POINT3D));
1426 
1427  p1.x = 1.0;
1428  p2.y = 1.0;
1429  angle = M_PI/4;
1430  vector_rotate(&p1, &p2, angle, &n);
1431  //printf("%g %g %g\n\n", n.x, n.y, n.z);
1432  CU_ASSERT_DOUBLE_EQUAL(n.x, 0.707107, 0.00001);
1433 
1434  angle = 2*M_PI/400000000;
1435  vector_rotate(&p1, &p2, angle, &n);
1436  //printf("%.21g %.21g %.21g\n\n", n.x, n.y, n.z);
1437  CU_ASSERT_DOUBLE_EQUAL(n.x, 0.999999999999999888978, 0.0000000000000001);
1438  CU_ASSERT_DOUBLE_EQUAL(n.y, 1.57079632679489654446e-08, 0.0000000000000001);
1439 
1440  angle = 0;
1441  vector_rotate(&p1, &p2, angle, &n);
1442  //printf("%.16g %.16g %.16g\n\n", n.x, n.y, n.z);
1443  CU_ASSERT_DOUBLE_EQUAL(n.x, 1.0, 0.00000001);
1444 }
1445 
1447 {
1448  LWGEOM *lwg1, *lwg2;
1449  LWLINE *lwl;
1450  double max = 100000.0 / WGS84_RADIUS;
1451  //char *wkt;
1452 
1453  /* Simple case */
1454  lwg1 = lwgeom_from_wkt("LINESTRING(0 20, 5 20)", LW_PARSER_CHECK_NONE);
1455  lwg2 = lwgeom_segmentize_sphere(lwg1, max);
1456  lwl = (LWLINE*)lwg2;
1457  //wkt = lwgeom_to_ewkt(lwg2);
1458  CU_ASSERT_EQUAL(lwl->points->npoints, 7);
1459  lwgeom_free(lwg1);
1460  lwgeom_free(lwg2);
1461  //lwfree(wkt);
1462 
1463  return;
1464 }
1465 
1466 static void test_lwgeom_area_sphere(void)
1467 {
1468  LWGEOM *lwg;
1469  double area;
1470  SPHEROID s;
1471 
1472  /* Init to WGS84 */
1474 
1475  /* Simple case */
1476  lwg = lwgeom_from_wkt("POLYGON((1 1, 1 2, 2 2, 2 1, 1 1))", LW_PARSER_CHECK_NONE);
1477  area = lwgeom_area_sphere(lwg, &s);
1478 
1479  CU_ASSERT_DOUBLE_EQUAL(area, 12360265021.3561, 1.0);
1480  lwgeom_free(lwg);
1481 
1482  /* Robustness tests, from ticket #3393 */
1483  lwg = lwgeom_from_wkt("POLYGON((0 78.703946026663,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026663))", LW_PARSER_CHECK_NONE);
1484  area = lwgeom_area_sphere(lwg, &s);
1485  CU_ASSERT_DOUBLE_EQUAL(area, 127516467322130, 1.0);
1486  lwgeom_free(lwg);
1487 
1488  lwg = lwgeom_from_wkt("POLYGON((0 78.703946026662,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026662))", LW_PARSER_CHECK_NONE);
1489  area = lwgeom_area_sphere(lwg, &s);
1490  CU_ASSERT_DOUBLE_EQUAL(area, 127516467322130, 1.0);
1491  lwgeom_free(lwg);
1492 
1493  lwg = lwgeom_from_wkt("POLYGON((0 78.703946026664,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026664))", LW_PARSER_CHECK_NONE);
1494  area = lwgeom_area_sphere(lwg, &s);
1495  CU_ASSERT_DOUBLE_EQUAL(area, 127516467322130, 1.0);
1496  lwgeom_free(lwg);
1497  /* end #3393 */
1498 }
1499 
1500 /*
1501 ** Used by test harness to register the tests in this file.
1502 */
1503 void geodetic_suite_setup(void);
1505 {
1506  CU_pSuite suite = CU_add_suite("Geodetic", NULL, NULL);
1510  PG_ADD_TEST(suite, test_signum);
1513  PG_ADD_TEST(suite, test_clairaut);
1524  PG_ADD_TEST(suite, test_gbox_utils);
1530 }
uint32_t gserialized_get_type(const GSERIALIZED *s)
Extract the geometry type from the serialized form (it hides in the anonymous data area...
Definition: g_serialized.c:56
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:2612
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:583
static void edge_deg2rad(GEOGRAPHIC_EDGE *e)
Convert an edge from degrees to radians.
Definition: cu_geodetic.c:27
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:897
static void test_ptarray_contains_point_sphere(void)
Definition: cu_geodetic.c:907
double longitude_radians_normalize(double lon)
Convert a longitude to the range of -PI,PI.
Definition: lwgeodetic.c:27
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:1261
void spheroid_init(SPHEROID *s, double a, double b)
Initialize a spheroid object for use in geodetic functions.
Definition: lwspheroid.c:20
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:1216
Two-point great circle segment from a to b.
Definition: lwgeodetic.h:42
#define PIR_B_TOUCH_RIGHT
Definition: lwgeodetic.h:77
void normalize(POINT3D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:564
static void test_sphere_direction(void)
Definition: cu_geodetic.c:71
void lwfree(void *mem)
Definition: lwutil.c:190
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
double y
Definition: liblwgeom.h:296
char * gbox_to_string(const GBOX *gbox)
Allocate a string representation of the GBOX, based on dimensionality of flags.
Definition: g_box.c:328
#define WGS84_RADIUS
Definition: liblwgeom.h:98
int npoints
Definition: liblwgeom.h:327
static void test_lwgeom_segmentize_sphere(void)
Definition: cu_geodetic.c:1446
#define PIR_A_TOUCH_RIGHT
Definition: lwgeodetic.h:75
#define POLYGONTYPE
Definition: liblwgeom.h:62
Datum area(PG_FUNCTION_ARGS)
double b
Definition: liblwgeom.h:270
uint8_t flags
Definition: liblwgeom.h:353
double xmax
Definition: liblwgeom.h:249
static void line2pts(const char *wkt, POINT3D *A1, POINT3D *A2)
Definition: cu_geodetic.c:541
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:2715
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1006
void lwline_free(LWLINE *line)
Definition: lwline.c:63
static void point_deg2rad(GEOGRAPHIC_POINT *p)
Convert an edge from radians to degrees.
Definition: cu_geodetic.c:49
#define signum(a)
Ape a java function.
Definition: lwgeodetic.h:66
double radius
Definition: liblwgeom.h:274
static void edge_set(double lon1, double lat1, double lon2, double lat2, GEOGRAPHIC_EDGE *e)
Definition: cu_geodetic.c:327
int gbox_centroid(const GBOX *gbox, POINT2D *out)
Computes the average(ish) center of the box and returns success.
Definition: lwgeodetic.c:244
double x
Definition: liblwgeom.h:296
static void test_lwgeom_area_sphere(void)
Definition: cu_geodetic.c:1466
#define WGS84_MINOR_AXIS
Definition: liblwgeom.h:97
static void test_sphere_project(void)
Definition: cu_geodetic.c:90
static void test_clairaut(void)
Definition: cu_geodetic.c:343
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition: lwin_wkt.c:844
#define FLAGS_SET_GEODETIC(flags, value)
Definition: liblwgeom.h:115
char gbox_data[][512]
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:33
static void test_edge_distance_to_edge(void)
Definition: cu_geodetic.c:741
static void test_lwgeom_check_geodetic(void)
Definition: cu_geodetic.c:780
char ** result
Definition: liblwgeom.h:218
static void test_ptarray_contains_point_sphere_iowa(void)
Definition: cu_geodetic.c:1122
static void test_gbox_utils(void)
Definition: cu_geodetic.c:1344
static void test_signum(void)
Definition: cu_geodetic.c:64
static void test_vector_angle(void)
Definition: cu_geodetic.c:1387
void geodetic_suite_setup(void)
Definition: cu_geodetic.c:1504
double z
Definition: liblwgeom.h:296
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate bounding box of a geometry, automatically taking into account whether it is cartesian or ge...
Definition: lwgeom.c:608
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:1706
double x
Definition: liblwgeom.h:284
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:1048
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:490
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:3201
double zmax
Definition: liblwgeom.h:253
double ymin
Definition: liblwgeom.h:250
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:249
double xmin
Definition: liblwgeom.h:248
#define LW_FALSE
Definition: liblwgeom.h:52
#define rad2deg(r)
Definition: lwgeodetic.h:61
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:44
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:51
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
Definition: lwgeodetic.c:1165
static void point_rad2deg(GEOGRAPHIC_POINT *p)
Convert a point from radians to degrees.
Definition: cu_geodetic.c:58
void ll2cart(const POINT2D *g, POINT3D *p)
Convert lon/lat coordinates to cartesion coordinates on unit sphere.
Definition: lwgeodetic.c:374
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:1023
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:154
#define PG_ADD_TEST(suite, testfunc)
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:522
POINTARRAY * ptarray_construct_reference_data(char hasz, char hasm, uint32_t npoints, uint8_t *ptlist)
Construct a new POINTARRAY, referencing to the data from ptlist.
Definition: ptarray.c:280
GSERIALIZED * gserialized_from_lwgeom(LWGEOM *geom, int is_geodetic, size_t *size)
Allocate a new GSERIALIZED from an LWGEOM.
Definition: g_serialized.c:908
static void test_lwpoly_covers_point2d(void)
Definition: cu_geodetic.c:1096
POINTARRAY ** rings
Definition: liblwgeom.h:413
#define RANDOM_TEST
Definition: cu_geodetic.c:22
static void test_gbox_from_spherical_coordinates(void)
Definition: cu_geodetic.c:204
double ymax
Definition: liblwgeom.h:251
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:45
double y
Definition: liblwgeom.h:284
char * s
Definition: cu_in_wkt.c:24
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:434
static void point_set(double lon, double lat, GEOGRAPHIC_POINT *p)
Definition: cu_geodetic.c:336
double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the geodetic area of a lwgeom on the sphere.
Definition: lwgeodetic.c:1924
uint8_t flags
Definition: liblwgeom.h:247
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:131
LWGEOM * lwgeom_from_hexwkb(const char *hexwkb, const char check)
Definition: lwin_wkb.c:753
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesion coordinates on unit sphere.
Definition: lwgeodetic.c:355
int gbox_data_length
int gbox_geocentric_slow
For testing geodetic bounding box, we have a magic global variable.
Definition: lwgeodetic.c:22
double latitude_radians_normalize(double lat)
Convert a latitude to the range of -PI/2,PI/2.
Definition: lwgeodetic.c:55
static void test_gserialized_from_lwgeom(void)
Definition: cu_geodetic.c:868
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:3096
static void test_spheroid_distance(void)
Definition: cu_geodetic.c:1240
static void test_edge_intersects(void)
Definition: cu_geodetic.c:563
double a
Definition: liblwgeom.h:269
#define WGS84_MAJOR_AXIS
Definition: liblwgeom.h:95
LWGEOM * lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length)
Derive a new geometry with vertices added to ensure no vertex is more than max_seg_length (in radians...
Definition: lwgeodetic.c:1638
static void test_edge_distance_to_point(void)
Definition: cu_geodetic.c:704
char * iowa_data
double zmin
Definition: liblwgeom.h:252
double gbox_angular_height(const GBOX *gbox)
GBOX utility functions to figure out coverage/location on the globe.
Definition: lwgeodetic.c:165
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:157
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:60
double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the geodetic area of a lwgeom on the spheroid.
Definition: lwspheroid.c:513
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:59
double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance)
Calculate the geodetic distance from lwgeom1 to lwgeom2 on the spheroid.
Definition: lwgeodetic.c:2078
static void test_edge_intersection(void)
Definition: cu_geodetic.c:381
static void test_gserialized_get_gbox_geocentric(void)
Definition: cu_geodetic.c:275
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:2384
#define PIR_B_TOUCH_LEFT
Definition: lwgeodetic.h:78
static LWGEOM * lwgeom_over_gserialized(char *wkt, GSERIALIZED **g)
Definition: cu_geodetic.c:769
static void test_vector_rotate(void)
Definition: cu_geodetic.c:1418
#define MULTILINETYPE
Definition: liblwgeom.h:64
double vector_angle(const POINT3D *v1, const POINT3D *v2)
Angle between two unit vectors.
Definition: lwgeodetic.c:454
#define PIR_COLINEAR
Definition: lwgeodetic.h:74
static void test_spheroid_area(void)
Definition: cu_geodetic.c:1281
static void test_lwgeom_distance_sphere(void)
Definition: cu_geodetic.c:1142
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:1074
#define PIR_INTERSECTS
Definition: lwgeodetic.h:73
POINTARRAY * points
Definition: liblwgeom.h:378
double gbox_angular_width(const GBOX *gbox)
Returns the angular width (longitudinal span) of the box in radians.
Definition: lwgeodetic.c:192
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:924