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