PostGIS  2.5.7dev-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 
32 static void edge_deg2rad(GEOGRAPHIC_EDGE *e)
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 
54 static void point_deg2rad(GEOGRAPHIC_POINT *p)
55 {
56  p->lat = latitude_radians_normalize(deg2rad(p->lat));
57  p->lon = longitude_radians_normalize(deg2rad(p->lon));
58 }
59 
63 static void point_rad2deg(GEOGRAPHIC_POINT *p)
64 {
65  p->lat = rad2deg(p->lat);
66  p->lon = rad2deg(p->lon);
67 }
68 
69 static void test_sphere_direction(void)
70 {
71  GEOGRAPHIC_POINT s, e;
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 {
94  GEOGRAPHIC_POINT s, e;
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
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 
224 static void test_gbox_from_spherical_coordinates(void)
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 
246  lwline = lwline_as_lwgeom(lwline_construct(SRID_UNKNOWN, 0, pa));
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 
263  gbox_geocentric_slow = LW_FALSE;
264  lwgeom_calculate_gbox_geodetic(lwline, &gbox);
265  gbox_geocentric_slow = LW_TRUE;
266  lwgeom_calculate_gbox_geodetic(lwline, &gbox_slow);
267  gbox_geocentric_slow = LW_FALSE;
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 
295 static void test_gserialized_get_gbox_geocentric(void)
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
308  lwg = lwgeom_from_wkt(gbox_data[i], LW_PARSER_CHECK_NONE);
309  FLAGS_SET_GEODETIC(lwg->flags, 1);
310  gbox_geocentric_slow = LW_FALSE;
311  lwgeom_calculate_gbox(lwg, &gbox);
312  gbox_geocentric_slow = LW_TRUE;
313  lwgeom_calculate_gbox(lwg, &gbox_slow);
314  gbox_geocentric_slow = LW_FALSE;
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;
404  GEOGRAPHIC_POINT g;
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 {
563  LWLINE *l = (LWLINE*)lwgeom_from_wkt(wkt, LW_PARSER_CHECK_NONE);
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;
586  GEOGRAPHIC_POINT g;
587  uint32_t rv;
588 
589  /* 5m close case */
590  line2pts("LINESTRING(58.5112113206308 0, 58.511211320077201 0.00090193752520337797)", &A1, &A2);
591  line2pts("LINESTRING(58.511166525601702 0.00027058124084120699, 58.511166525562899 0.00036077498778824899)", &B1, &B2);
592  rv = edge_intersects(&A1, &A2, &B1, &B2);
593  CU_ASSERT(rv == 0);
594 
595  /* Covers case, end-to-end intersection */
596  line2pts("LINESTRING(50 -10.999999999999998224, -10.0 50.0)", &A1, &A2);
597  line2pts("LINESTRING(-10.0 50.0, -10.272779983831613393 -16.937003313332997578)", &B1, &B2);
598  rv = edge_intersects(&A1, &A2, &B1, &B2);
599  CU_ASSERT(rv & PIR_INTERSECTS);
600 
601  /* Medford case, very short segment vs very long one */
602  g.lat = 0.74123572595649878103;
603  g.lon = -2.1496353191142714145;
604  geog2cart(&g, &A1);
605  g.lat = 0.74123631950116664058;
606  g.lon = -2.1496353248304860273;
607  geog2cart(&g, &A2);
608  g.lat = 0.73856343764436815924;
609  g.lon = -2.1461493501950630325;
610  geog2cart(&g, &B1);
611  g.lat = 0.70971354024834598651;
612  g.lon = 2.1082194552519770703;
613  geog2cart(&g, &B2);
614  rv = edge_intersects(&A1, &A2, &B1, &B2);
615  CU_ASSERT(rv == 0);
616 
617  /* Second Medford case, very short segment vs very long one */
618  g.lat = 0.73826546728290887156;
619  g.lon = -2.14426380171833042;
620  geog2cart(&g, &A1);
621  g.lat = 0.73826545883786642843;
622  g.lon = -2.1442638997530165668;
623  geog2cart(&g, &A2);
624  g.lat = 0.73775469118192538165;
625  g.lon = -2.1436035534281718817;
626  geog2cart(&g, &B1);
627  g.lat = 0.71021099548296817705;
628  g.lon = 2.1065275171200439353;
629  geog2cart(&g, &B2);
630  rv = edge_intersects(&A1, &A2, &B1, &B2);
631  CU_ASSERT(rv == PIR_INTERSECTS);
632 
633  /* Again, this time with a less exact input edge. */
634  line2pts("LINESTRING(-123.165031277506 42.4696787216231, -123.165031605021 42.4697127292275)", &A1, &A2);
635  rv = edge_intersects(&A1, &A2, &B1, &B2);
636  CU_ASSERT(rv == 0);
637 
638  /* 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 1.0)", &B1, &B2);
641  rv = edge_intersects(&A1, &A2, &B1, &B2);
642  CU_ASSERT(rv == PIR_INTERSECTS);
643 
644  /* No intersection 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 -2.0)", &B1, &B2);
647  rv = edge_intersects(&A1, &A2, &B1, &B2);
648  CU_ASSERT(rv == 0);
649 
650  /* End touches middle of segment at (0 0) */
651  line2pts("LINESTRING(-1.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) );
655 
656  /* End touches end of segment at (0 0) */
657  line2pts("LINESTRING(0.0 0.0, 1.0 0.0)", &A1, &A2);
658  line2pts("LINESTRING(0.0 -1.0, 0.0 0.0)", &B1, &B2);
659  rv = edge_intersects(&A1, &A2, &B1, &B2);
660  CU_ASSERT(rv == (PIR_INTERSECTS|PIR_B_TOUCH_RIGHT|PIR_A_TOUCH_RIGHT) );
661 
662  /* Intersection at (180 0) */
663  line2pts("LINESTRING(-179.0 -1.0, 179.0 1.0)", &A1, &A2);
664  line2pts("LINESTRING(-179.0 1.0, 179.0 -1.0)", &B1, &B2);
665  rv = edge_intersects(&A1, &A2, &B1, &B2);
666  CU_ASSERT(rv == PIR_INTERSECTS);
667 
668  /* Intersection at (180 0) */
669  line2pts("LINESTRING(-170.0 0.0, 170.0 0.0)", &A1, &A2);
670  line2pts("LINESTRING(180.0 -10.0, 180.0 10.0)", &B1, &B2);
671  rv = edge_intersects(&A1, &A2, &B1, &B2);
672  CU_ASSERT(rv == PIR_INTERSECTS);
673 
674  /* Intersection at north pole */
675  line2pts("LINESTRING(-180.0 80.0, 0.0 80.0)", &A1, &A2);
676  line2pts("LINESTRING(90.0 80.0, -90.0 80.0)", &B1, &B2);
677  rv = edge_intersects(&A1, &A2, &B1, &B2);
678  CU_ASSERT(rv == PIR_INTERSECTS);
679 
680  /* Equal edges return true */
681  line2pts("LINESTRING(45.0 10.0, 50.0 20.0)", &A1, &A2);
682  line2pts("LINESTRING(45.0 10.0, 50.0 20.0)", &B1, &B2);
683  rv = edge_intersects(&A1, &A2, &B1, &B2);
684  CU_ASSERT(rv & PIR_INTERSECTS);
685 
686  /* Parallel edges (same great circle, different end points) return true */
687  line2pts("LINESTRING(40.0 0.0, 70.0 0.0)", &A1, &A2);
688  line2pts("LINESTRING(60.0 0.0, 50.0 0.0)", &B1, &B2);
689  rv = edge_intersects(&A1, &A2, &B1, &B2);
690  CU_ASSERT(rv == (PIR_INTERSECTS|PIR_COLINEAR) );
691 
692  /* End touches arc at north pole */
693  line2pts("LINESTRING(-180.0 80.0, 0.0 80.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) );
697 
698  /* End touches end at north pole */
699  line2pts("LINESTRING(-180.0 80.0, 0.0 90.0)", &A1, &A2);
700  line2pts("LINESTRING(90.0 80.0, -90.0 90.0)", &B1, &B2);
701  rv = edge_intersects(&A1, &A2, &B1, &B2);
702  CU_ASSERT(rv == (PIR_INTERSECTS|PIR_B_TOUCH_LEFT|PIR_A_TOUCH_RIGHT) );
703 
704  /* Antipodal straddles. Great circles cross but at opposite */
705  /* sides of the globe */
706  /* #2534 */
707  /* http://www.gcmap.com/mapui?P=60N+90E-20S+90E%0D%0A0N+0E-90.04868865037885W+57.44011727050777S%0D%0A&MS=wls&DU=mi */
708  line2pts("LINESTRING(90.0 60.0, 90.0 -20.0)", &A1, &A2);
709  line2pts("LINESTRING(0.0 0.0, -90.04868865037885 -57.44011727050777)", &B1, &B2);
710  rv = edge_intersects(&A1, &A2, &B1, &B2);
711  CU_ASSERT(rv == 0);
712 
713  line2pts("LINESTRING(-5 0, 5 0)", &A1, &A2);
714  line2pts("LINESTRING(179 -5, 179 5)", &B1, &B2);
715  rv = edge_intersects(&A1, &A2, &B1, &B2);
716  CU_ASSERT(rv == 0);
717 
718  line2pts("LINESTRING(175 -85, 175 85)", &A1, &A2);
719  line2pts("LINESTRING(65 0, -105 0)", &B1, &B2);
720  rv = edge_intersects(&A1, &A2, &B1, &B2);
721  CU_ASSERT(rv == 0);
722 
723  line2pts("LINESTRING(175 -85, 175 85)", &A1, &A2);
724  line2pts("LINESTRING(45 0, -125 0)", &B1, &B2);
725  rv = edge_intersects(&A1, &A2, &B1, &B2);
726  CU_ASSERT(rv == 0);
727 
728 }
729 
730 static void test_edge_distance_to_point(void)
731 {
732  GEOGRAPHIC_EDGE e;
733  GEOGRAPHIC_POINT g;
734  GEOGRAPHIC_POINT closest;
735  double d;
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, 1.0, &g);
740  d = edge_distance_to_point(&e, &g, 0);
741  CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001);
742 
743  /* closest point at origin, one degree away */
744  edge_set(-50.0, 0.0, 50.0, 0.0, &e);
745  point_set(0.0, 2.0, &g);
746  d = edge_distance_to_point(&e, &g, &closest);
747 #if 0
748  printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e.start.lon, e.start.lat, e.end.lon, e.end.lat);
749  printf("POINT(%.9g %.9g)\n", g.lon, g.lat);
750  printf("\nDISTANCE == %.8g\n", d);
751 #endif
752  CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 90.0, 0.00001);
753  CU_ASSERT_DOUBLE_EQUAL(closest.lat, 0.0, 0.00001);
754  CU_ASSERT_DOUBLE_EQUAL(closest.lon, 0.0, 0.00001);
755 
756  /* Ticket #2351 */
757  edge_set(149.386990599235, -26.3567415843982, 149.386990599247, -26.3567415843965, &e);
758  point_set(149.386990599235, -26.3567415843982, &g);
759  d = edge_distance_to_point(&e, &g, &closest);
760  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
761  // printf("CLOSE POINT(%g %g)\n", closest.lon, closest.lat);
762  // printf(" ORIG POINT(%g %g)\n", g.lon, g.lat);
763  CU_ASSERT_DOUBLE_EQUAL(g.lat, closest.lat, 0.00001);
764  CU_ASSERT_DOUBLE_EQUAL(g.lon, closest.lon, 0.00001);
765 }
766 
767 static void test_edge_distance_to_edge(void)
768 {
769  GEOGRAPHIC_EDGE e1, e2;
770  GEOGRAPHIC_POINT c1, c2;
771  double d;
772 
773  /* closest point at origin, one degree away */
774  edge_set(-50.0, 0.0, 50.0, 0.0, &e1);
775  edge_set(-5.0, 20.0, 0.0, 1.0, &e2);
776  d = edge_distance_to_edge(&e1, &e2, &c1, &c2);
777 #if 0
778  printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e1.start.lon, e1.start.lat, e1.end.lon, e1.end.lat);
779  printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e2.start.lon, e2.start.lat, e2.end.lon, e2.end.lat);
780  printf("\nDISTANCE == %.8g\n", d);
781 #endif
782  CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001);
783  CU_ASSERT_DOUBLE_EQUAL(c1.lat, 0.0, 0.00001);
784  CU_ASSERT_DOUBLE_EQUAL(c2.lat, M_PI / 180.0, 0.00001);
785  CU_ASSERT_DOUBLE_EQUAL(c1.lon, 0.0, 0.00001);
786  CU_ASSERT_DOUBLE_EQUAL(c2.lon, 0.0, 0.00001);
787 }
788 
789 
790 
791 /*
792 * Build LWGEOM on top of *aligned* structure so we can use the read-only
793 * point access methods on them.
794 */
795 static LWGEOM* lwgeom_over_gserialized(char *wkt, GSERIALIZED **g)
796 {
797  LWGEOM *lwg;
798 
799  lwg = lwgeom_from_wkt(wkt, LW_PARSER_CHECK_NONE);
800  FLAGS_SET_GEODETIC(lwg->flags, 1);
801  *g = gserialized_from_lwgeom(lwg, 0);
802  lwgeom_free(lwg);
803  return lwgeom_from_gserialized(*g);
804 }
805 
806 static void test_lwgeom_check_geodetic(void)
807 {
808  LWGEOM *geom;
809  int i = 0;
810 
811  char ewkt[][512] =
812  {
813  "POINT(0 0.2)",
814  "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
815  "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))",
816  "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
817  "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)))",
818  "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))))",
819  "POINT(0 220.2)",
820  "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
821  "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))",
822  "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
823  "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)))",
824  "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))))",
825  };
826 
827  for ( i = 0; i < 6; i++ )
828  {
829  GSERIALIZED *g;
830  geom = lwgeom_over_gserialized(ewkt[i], &g);
831  CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), LW_TRUE);
832  lwgeom_free(geom);
833  lwfree(g);
834  }
835 
836  for ( i = 6; i < 12; i++ )
837  {
838  GSERIALIZED *g;
839  //char *out_ewkt;
840  geom = lwgeom_over_gserialized(ewkt[i], &g);
841  CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), LW_FALSE);
842  //out_ewkt = lwgeom_to_ewkt(geom);
843  //printf("%s\n", out_ewkt);
844  lwgeom_free(geom);
845  lwfree(g);
846  }
847 
848 }
849 
850 /*
851 static void test_gbox_calculation(void)
852 {
853 
854  LWGEOM *geom;
855  int i = 0;
856  GBOX *gbox = gbox_new(gflags(0,0,0));
857  BOX3D *box3d;
858 
859  char ewkt[][512] =
860  {
861  "POINT(0 0.2)",
862  "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
863  "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))",
864  "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
865  "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)))",
866  "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))))",
867  "POINT(0 220.2)",
868  "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
869  "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))",
870  "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
871  "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)))",
872  "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))))",
873  };
874 
875  for ( i = 0; i < 6; i++ )
876  {
877  GSERIALIZED *g;
878  geom = lwgeom_over_gserialized(ewkt[i], &g);
879  lwgeom_calculate_gbox_cartesian(geom, gbox);
880  box3d = lwgeom_compute_box3d(geom);
881  //printf("%g %g\n", gbox->xmin, box3d->xmin);
882  CU_ASSERT_EQUAL(gbox->xmin, box3d->xmin);
883  CU_ASSERT_EQUAL(gbox->xmax, box3d->xmax);
884  CU_ASSERT_EQUAL(gbox->ymin, box3d->ymin);
885  CU_ASSERT_EQUAL(gbox->ymax, box3d->ymax);
886  lwgeom_free(geom);
887  lwfree(box3d);
888  lwfree(g);
889  }
890  lwfree(gbox);
891 }
892 */
893 
894 static void test_gserialized_from_lwgeom(void)
895 {
896  LWGEOM *geom;
897  GSERIALIZED *g;
898  uint32_t type;
899  double *inspect; /* To poke right into the blob. */
900 
901  geom = lwgeom_from_wkt("POINT(0 0.2)", LW_PARSER_CHECK_NONE);
902  FLAGS_SET_GEODETIC(geom->flags, 1);
903  g = gserialized_from_lwgeom(geom, 0);
904  type = gserialized_get_type(g);
905  CU_ASSERT_EQUAL( type, POINTTYPE );
906  inspect = (double*)g;
907  CU_ASSERT_EQUAL(inspect[3], 0.2);
908  lwgeom_free(geom);
909  lwfree(g);
910 
911  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);
912  FLAGS_SET_GEODETIC(geom->flags, 1);
913  g = gserialized_from_lwgeom(geom, 0);
914  type = gserialized_get_type(g);
915  CU_ASSERT_EQUAL( type, POLYGONTYPE );
916  inspect = (double*)g;
917  CU_ASSERT_EQUAL(inspect[9], 2.5);
918  lwgeom_free(geom);
919  lwfree(g);
920 
921  geom = lwgeom_from_wkt("MULTILINESTRING((0 0, 1 1),(0 0.1, 1 1))", LW_PARSER_CHECK_NONE);
922  FLAGS_SET_GEODETIC(geom->flags, 1);
923  g = gserialized_from_lwgeom(geom, 0);
924  type = gserialized_get_type(g);
925  CU_ASSERT_EQUAL( type, MULTILINETYPE );
926  inspect = (double*)g;
927  CU_ASSERT_EQUAL(inspect[12], 0.1);
928  lwgeom_free(geom);
929  lwfree(g);
930 
931 }
932 
933 static void test_ptarray_contains_point_sphere(void)
934 {
935  LWGEOM *lwg;
936  LWPOLY *poly;
937  POINT2D pt_to_test;
938  POINT2D pt_outside;
939  int result;
940 
941  /* Small polygon and huge distance between outside point and close-but-not-quite-inside point. Should return LW_FALSE. Pretty degenerate case. */
942  lwg = lwgeom_from_hexwkb("0103000020E61000000100000025000000ACAD6F91DDB65EC03F84A86D57264540CCABC279DDB65EC0FCE6926B57264540B6DEAA62DDB65EC0A79F6B63572645402E0BE84CDDB65EC065677155572645405D0B1D39DDB65EC0316310425726454082B5DB27DDB65EC060A4E12957264540798BB619DDB65EC0C393A10D57264540D4BC160FDDB65EC0BD0320EE56264540D7AC4E08DDB65EC096C862CC56264540AFD29205DDB65EC02A1F68A956264540363AFA06DDB65EC0722E418656264540B63A780CDDB65EC06E9B0064562645409614E215DDB65EC0E09DA84356264540FF71EF22DDB65EC0B48145265626454036033F33DDB65EC081B8A60C5626454066FB4546DDB65EC08A47A6F7552645409061785BDDB65EC0F05AE0E755264540D4B63772DDB65EC05C86CEDD55264540D2E4C689DDB65EC09B6EBFD95526454082E573A1DDB65EC0C90BD5DB552645401ABE85B8DDB65EC06692FCE35526454039844ECEDDB65EC04D8AF6F155264540928319E2DDB65EC0AD8D570556264540D31055F3DDB65EC02D618F1D56264540343B7A01DEB65EC0EB70CF3956264540920A1A0CDEB65EC03B00515956264540911BE212DEB65EC0E43A0E7B56264540E3F69D15DEB65EC017E4089E562645408D903614DEB65EC0F0D42FC1562645402191B80EDEB65EC0586870E35626454012B84E05DEB65EC09166C80357264540215B41F8DDB65EC08F832B21572645408392F7E7DDB65EC01138C13A57264540F999F0D4DDB65EC0E4A9C14F57264540AC3FB8BFDDB65EC0EED6875F57264540D3DCFEA8DDB65EC04F6C996957264540ACAD6F91DDB65EC03F84A86D57264540", LW_PARSER_CHECK_NONE);
943  poly = (LWPOLY*)lwg;
944  pt_to_test.x = -122.819436560680316;
945  pt_to_test.y = 42.2702301207017328;
946  pt_outside.x = 120.695136159150778;
947  pt_outside.y = 40.6920926049588516;
948  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
949  CU_ASSERT_EQUAL(result, LW_FALSE);
950  lwgeom_free(lwg);
951 
952  /* Point on ring between vertexes case */
953  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);
954  poly = (LWPOLY*)lwg;
955  pt_to_test.x = 1.1;
956  pt_to_test.y = 1.05;
957  pt_outside.x = 1.2;
958  pt_outside.y = 1.05;
959  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
960  CU_ASSERT_EQUAL(result, LW_TRUE);
961  lwgeom_free(lwg);
962 
963  /* Simple containment case */
964  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);
965  poly = (LWPOLY*)lwg;
966  pt_to_test.x = 1.05;
967  pt_to_test.y = 1.05;
968  pt_outside.x = 1.2;
969  pt_outside.y = 1.15;
970  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
971  CU_ASSERT_EQUAL(result, LW_TRUE);
972  lwgeom_free(lwg);
973 
974  /* Less Simple containment case. */
975  /* Interior point quite close to boundary and stab line going through bottom edge vertex */
976  /* This breaks the "extend-it" trick of handling vertex crossings */
977  /* It should also break the "lowest end" trick. */
978  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);
979  poly = (LWPOLY*)lwg;
980  pt_to_test.x = 1.05;
981  pt_to_test.y = 1.00;
982  pt_outside.x = 1.05;
983  pt_outside.y = 0.5;
984  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
985  CU_ASSERT_EQUAL(result, LW_TRUE);
986  lwgeom_free(lwg);
987 
988  /* Simple noncontainment case */
989  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);
990  poly = (LWPOLY*)lwg;
991  pt_to_test.x = 1.05;
992  pt_to_test.y = 1.15;
993  pt_outside.x = 1.2;
994  pt_outside.y = 1.2;
995  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
996  CU_ASSERT_EQUAL(result, LW_FALSE);
997  lwgeom_free(lwg);
998 
999  /* Harder noncontainment case */
1000  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);
1001  poly = (LWPOLY*)lwg;
1002  pt_to_test.x = 1.05;
1003  pt_to_test.y = 0.9;
1004  pt_outside.x = 1.2;
1005  pt_outside.y = 1.05;
1006  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1007  CU_ASSERT_EQUAL(result, LW_FALSE);
1008  lwgeom_free(lwg);
1009 
1010  /* Harder containment case */
1011  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);
1012  poly = (LWPOLY*)lwg;
1013  pt_to_test.x = 1.0;
1014  pt_to_test.y = 1.0;
1015  pt_outside.x = 1.0;
1016  pt_outside.y = 10.0;
1017  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1018  CU_ASSERT_EQUAL(result, LW_TRUE);
1019  lwgeom_free(lwg);
1020 
1021  /* Point on ring at first vertex case */
1022  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);
1023  poly = (LWPOLY*)lwg;
1024  pt_to_test.x = 1.0;
1025  pt_to_test.y = 1.0;
1026  pt_outside.x = 1.2;
1027  pt_outside.y = 1.05;
1028  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1029  CU_ASSERT_EQUAL(result, LW_TRUE);
1030  lwgeom_free(lwg);
1031 
1032  /* Point on ring at vertex case */
1033  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);
1034  poly = (LWPOLY*)lwg;
1035  pt_to_test.x = 1.0;
1036  pt_to_test.y = 1.1;
1037  pt_outside.x = 1.2;
1038  pt_outside.y = 1.05;
1039  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1040  CU_ASSERT_EQUAL(result, LW_TRUE);
1041  lwgeom_free(lwg);
1042 
1043  /* Co-linear crossing case for point-in-polygon test, should return LW_TRUE */
1044  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);
1045  poly = (LWPOLY*)lwg;
1046  pt_to_test.x = 1.1;
1047  pt_to_test.y = 1.05;
1048  pt_outside.x = 1.1;
1049  pt_outside.y = 1.3;
1050  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1051  CU_ASSERT_EQUAL(result, LW_TRUE);
1052  lwgeom_free(lwg);
1053 
1054  /* Co-linear grazing case for point-in-polygon test, should return LW_FALSE */
1055  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);
1056  poly = (LWPOLY*)lwg;
1057  pt_to_test.x = 1.0;
1058  pt_to_test.y = 0.0;
1059  pt_outside.x = 1.0;
1060  pt_outside.y = 2.0;
1061  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1062  CU_ASSERT_EQUAL(result, LW_FALSE);
1063  lwgeom_free(lwg);
1064 
1065  /* Grazing case for point-in-polygon test, should return LW_FALSE */
1066  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);
1067  poly = (LWPOLY*)lwg;
1068  pt_to_test.x = 1.5;
1069  pt_to_test.y = 1.0;
1070  pt_outside.x = 1.5;
1071  pt_outside.y = 2.0;
1072  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1073  CU_ASSERT_EQUAL(result, LW_FALSE);
1074  lwgeom_free(lwg);
1075 
1076  /* Grazing case at first point for point-in-polygon test, should return LW_FALSE */
1077  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);
1078  poly = (LWPOLY*)lwg;
1079  pt_to_test.x = 1.0;
1080  pt_to_test.y = 0.0;
1081  pt_outside.x = 1.0;
1082  pt_outside.y = 2.0;
1083  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1084  CU_ASSERT_EQUAL(result, LW_FALSE);
1085  lwgeom_free(lwg);
1086 
1087  /* Outside multi-crossing case for point-in-polygon test, should return LW_FALSE */
1088  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);
1089  poly = (LWPOLY*)lwg;
1090  pt_to_test.x = 0.99;
1091  pt_to_test.y = 0.99;
1092  pt_outside.x = 1.21;
1093  pt_outside.y = 1.21;
1094  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1095  CU_ASSERT_EQUAL(result, LW_FALSE);
1096  lwgeom_free(lwg);
1097 
1098  /* Inside multi-crossing case for point-in-polygon test, should return LW_TRUE */
1099  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);
1100  poly = (LWPOLY*)lwg;
1101  pt_to_test.x = 1.11;
1102  pt_to_test.y = 1.11;
1103  pt_outside.x = 1.21;
1104  pt_outside.y = 1.21;
1105  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1106  CU_ASSERT_EQUAL(result, LW_TRUE);
1107  lwgeom_free(lwg);
1108 
1109  /* Point on vertex of ring */
1110  lwg = lwgeom_from_wkt("POLYGON((-9 50,51 -11,-10 50,-9 50))", LW_PARSER_CHECK_NONE);
1111  poly = (LWPOLY*)lwg;
1112  pt_to_test.x = -10.0;
1113  pt_to_test.y = 50.0;
1114  pt_outside.x = -10.2727799838316134;
1115  pt_outside.y = -16.9370033133329976;
1116  result = ptarray_contains_point_sphere(poly->rings[0], &pt_outside, &pt_to_test);
1117  CU_ASSERT_EQUAL(result, LW_TRUE);
1118  lwgeom_free(lwg);
1119 
1120 }
1121 
1122 static void test_lwpoly_covers_point2d(void)
1123 {
1124  LWPOLY *poly;
1125  LWGEOM *lwg;
1126  POINT2D pt_to_test;
1127  int result;
1128 
1129  lwg = lwgeom_from_wkt("POLYGON((-9 50,51 -11,-10 50,-9 50))", LW_PARSER_CHECK_NONE);
1130  poly = (LWPOLY*)lwg;
1131  pt_to_test.x = -10.0;
1132  pt_to_test.y = 50.0;
1133  result = lwpoly_covers_point2d(poly, &pt_to_test);
1134  CU_ASSERT_EQUAL(result, LW_TRUE);
1135  lwgeom_free(lwg);
1136 
1137  /* Great big ring */
1138  lwg = lwgeom_from_wkt("POLYGON((-40.0 52.0, -67.0 -29.0, 102.0 -6.0, -40.0 52.0))", LW_PARSER_CHECK_NONE);
1139  poly = (LWPOLY*)lwg;
1140  pt_to_test.x = 4.0;
1141  pt_to_test.y = 11.0;
1142  result = lwpoly_covers_point2d(poly, &pt_to_test);
1143  CU_ASSERT_EQUAL(result, LW_TRUE);
1144  lwgeom_free(lwg);
1145 
1146  /* Triangle over the antimeridian */
1147  lwg = lwgeom_from_wkt("POLYGON((140 52, 152.0 -6.0, -120.0 -29.0, 140 52))", LW_PARSER_CHECK_NONE);
1148  poly = (LWPOLY*)lwg;
1149  pt_to_test.x = -172.0;
1150  pt_to_test.y = -13.0;
1151  result = lwpoly_covers_point2d(poly, &pt_to_test);
1152  CU_ASSERT_EQUAL(result, LW_TRUE);
1153  lwgeom_free(lwg);
1154 
1155 }
1156 
1157 static void test_ptarray_contains_point_sphere_iowa(void)
1158 {
1159  LWGEOM *lwg = lwgeom_from_wkt(iowa_data, LW_PARSER_CHECK_NONE);
1160  LWPOLY *poly = (LWPOLY*)lwg;
1161  POINTARRAY *pa = poly->rings[0];
1162  POINT2D pt_outside, pt_to_test;
1163  int rv;
1164 
1165  pt_to_test.x = -95.900000000000006;
1166  pt_to_test.y = 42.899999999999999;
1167  pt_outside.x = -96.381873780830645;
1168  pt_outside.y = 40.185394449416371;
1169 
1170  rv = ptarray_contains_point_sphere(pa, &pt_outside, &pt_to_test);
1171  CU_ASSERT_EQUAL(rv, LW_TRUE);
1172 
1173  lwgeom_free(lwg);
1174 }
1175 
1176 
1177 static void test_lwgeom_distance_sphere(void)
1178 {
1179  LWGEOM *lwg1, *lwg2;
1180  double d;
1181  SPHEROID s;
1182 
1183  /* Init and force spherical */
1184  spheroid_init(&s, 6378137.0, 6356752.314245179498);
1185  s.a = s.b = s.radius;
1186 
1187  /* Line/line distance, 1 degree apart */
1188  lwg1 = lwgeom_from_wkt("LINESTRING(-30 10, -20 5, -10 3, 0 1)", LW_PARSER_CHECK_NONE);
1189  lwg2 = lwgeom_from_wkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", LW_PARSER_CHECK_NONE);
1190  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1191  CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001);
1192  lwgeom_free(lwg1);
1193  lwgeom_free(lwg2);
1194 
1195  /* Line/line distance, crossing, 0.0 apart */
1196  lwg1 = lwgeom_from_wkt("LINESTRING(-30 10, -20 5, -10 3, 0 1)", LW_PARSER_CHECK_NONE);
1197  lwg2 = lwgeom_from_wkt("LINESTRING(-10 -5, -5 20, 5 0, 10 -5)", LW_PARSER_CHECK_NONE);
1198  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1199  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1200  lwgeom_free(lwg1);
1201  lwgeom_free(lwg2);
1202 
1203  /* Line/point distance, 1 degree apart */
1204  lwg1 = lwgeom_from_wkt("POINT(-4 1)", LW_PARSER_CHECK_NONE);
1205  lwg2 = lwgeom_from_wkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", LW_PARSER_CHECK_NONE);
1206  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1207  CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001);
1208  lwgeom_free(lwg1);
1209  lwgeom_free(lwg2);
1210 
1211  lwg1 = lwgeom_from_wkt("POINT(-4 1)", LW_PARSER_CHECK_NONE);
1212  lwg2 = lwgeom_from_wkt("POINT(-4 -1)", LW_PARSER_CHECK_NONE);
1213  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1214  CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 90.0, 0.00001);
1215  lwgeom_free(lwg1);
1216  lwgeom_free(lwg2);
1217 
1218  /* Poly/point distance, point inside polygon, 0.0 apart */
1219  lwg1 = lwgeom_from_wkt("POLYGON((-4 1, -3 5, 1 2, 1.5 -5, -4 1))", LW_PARSER_CHECK_NONE);
1220  lwg2 = lwgeom_from_wkt("POINT(-1 -1)", LW_PARSER_CHECK_NONE);
1221  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1222  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1223  lwgeom_free(lwg1);
1224  lwgeom_free(lwg2);
1225 
1226  /* Poly/point distance, point inside polygon hole, 1 degree apart */
1227  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);
1228  lwg2 = lwgeom_from_wkt("POINT(-1 -1)", LW_PARSER_CHECK_NONE);
1229  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1230  CU_ASSERT_DOUBLE_EQUAL(d, 111178.142466, 0.1);
1231  lwgeom_free(lwg1);
1232  lwgeom_free(lwg2);
1233 
1234  /* Poly/point distance, point on hole boundary, 0.0 apart */
1235  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);
1236  lwg2 = lwgeom_from_wkt("POINT(2 2)", LW_PARSER_CHECK_NONE);
1237  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1238  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1239  lwgeom_free(lwg1);
1240  lwgeom_free(lwg2);
1241 
1242  /* Medford test case #1 */
1243  lwg1 = lwgeom_from_hexwkb("0105000020E610000001000000010200000002000000EF7B8779C7BD5EC0FD20D94B852845400E539C62B9BD5EC0F0A5BE767C284540", LW_PARSER_CHECK_NONE);
1244  lwg2 = lwgeom_from_hexwkb("0106000020E61000000100000001030000000100000007000000280EC3FB8CCA5EC0A5CDC747233C45402787C8F58CCA5EC0659EA2761E3C45400CED58DF8FCA5EC0C37FAE6E1E3C4540AE97B8E08FCA5EC00346F58B1F3C4540250359FD8ECA5EC05460628E1F3C45403738F4018FCA5EC05DC84042233C4540280EC3FB8CCA5EC0A5CDC747233C4540", LW_PARSER_CHECK_NONE);
1245  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1246  CU_ASSERT_DOUBLE_EQUAL(d, 23630.8003, 0.1);
1247  lwgeom_free(lwg1);
1248  lwgeom_free(lwg2);
1249 
1250  /* Ticket #2351 */
1251  lwg1 = lwgeom_from_wkt("LINESTRING(149.386990599235 -26.3567415843982,149.386990599247 -26.3567415843965)", LW_PARSER_CHECK_NONE);
1252  lwg2 = lwgeom_from_wkt("POINT(149.386990599235 -26.3567415843982)", LW_PARSER_CHECK_NONE);
1253  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1254  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1255  lwgeom_free(lwg1);
1256  lwgeom_free(lwg2);
1257 
1258  /* Ticket #2638, no "M" */
1259  lwg1 = lwgeom_from_wkt("LINESTRING (-41.0821 50.3036,50 -41)", LW_PARSER_CHECK_NONE);
1260  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);
1261  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1262  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1263  lwgeom_free(lwg1);
1264  lwgeom_free(lwg2);
1265 
1266  /* Ticket #2638, with "M" */
1267  lwg1 = lwgeom_from_wkt("LINESTRING M (-41.0821 50.3036 1,50 -41 1)", LW_PARSER_CHECK_NONE);
1268  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);
1269  d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
1270  CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
1271  lwgeom_free(lwg1);
1272  lwgeom_free(lwg2);
1273 }
1274 
1275 static void test_spheroid_distance(void)
1276 {
1277  GEOGRAPHIC_POINT g1, g2;
1278  double d;
1279 #if ! PROJ_GEODESIC
1280  double epsilon; /* irregular */
1281 #else
1282  const double epsilon = 1e-8; /* at least 10 nm precision */
1283 #endif
1284  SPHEROID s;
1285 
1286  /* Init to WGS84 */
1287  spheroid_init(&s, 6378137.0, 6356752.314245179498);
1288 
1289  /* One vertical degree
1290  $ GeodSolve -E -i -p 16 --input-string "0 0 1 0" */
1291  point_set(0.0, 0.0, &g1);
1292  point_set(0.0, 1.0, &g2);
1293  d = spheroid_distance(&g1, &g2, &s);
1294 #if ! PROJ_GEODESIC
1295  epsilon = 1e-6;
1296 #endif
1297  CU_ASSERT_DOUBLE_EQUAL(d, 110574.3885577987957342, epsilon);
1298 
1299  /* Ten horizontal degrees
1300  $ GeodSolve -E -i -p 16 --input-string "0 -10 0 0" */
1301  point_set(-10.0, 0.0, &g1);
1302  point_set(0.0, 0.0, &g2);
1303  d = spheroid_distance(&g1, &g2, &s);
1304 #if ! PROJ_GEODESIC
1305  epsilon = 1e-3;
1306 #endif
1307  CU_ASSERT_DOUBLE_EQUAL(d, 1113194.9079327357264771, epsilon);
1308 
1309  /* One horizonal degree
1310  $ GeodSolve -E -i -p 16 --input-string "0 -1 0 0" */
1311  point_set(-1.0, 0.0, &g1);
1312  point_set(0.0, 0.0, &g2);
1313  d = spheroid_distance(&g1, &g2, &s);
1314 #if ! PROJ_GEODESIC
1315  epsilon = 1e-4;
1316 #endif
1317  CU_ASSERT_DOUBLE_EQUAL(d, 111319.4907932735726477, epsilon);
1318 
1319  /* Around world w/ slight bend
1320  $ GeodSolve -E -i -p 16 --input-string "0 -180 1 0" */
1321  point_set(-180.0, 0.0, &g1);
1322  point_set(0.0, 1.0, &g2);
1323  d = spheroid_distance(&g1, &g2, &s);
1324 #if ! PROJ_GEODESIC
1325  epsilon = 1e-5;
1326 #endif
1327  CU_ASSERT_DOUBLE_EQUAL(d, 19893357.0700676468277450, epsilon);
1328 
1329  /* Up to pole
1330  $ GeodSolve -E -i -p 16 --input-string "0 -180 90 0" */
1331  point_set(-180.0, 0.0, &g1);
1332  point_set(0.0, 90.0, &g2);
1333  d = spheroid_distance(&g1, &g2, &s);
1334 #if ! PROJ_GEODESIC
1335  epsilon = 1e-6;
1336 #endif
1337  CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7293127228117396, epsilon);
1338 
1339 }
1340 
1341 static void test_spheroid_area(void)
1342 {
1343  LWGEOM *lwg;
1344  GBOX gbox;
1345  double a1, a2;
1346  SPHEROID s;
1347 
1348  /* Init to WGS84 */
1349  spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
1350 
1351  gbox.flags = gflags(0, 0, 1);
1352 
1353  /* Medford lot test polygon */
1354  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);
1355  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1356  /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE -r --input-string \
1357  "42.5007249610493 -122.848227067007;42.5007179884263 -122.848309475585;"\
1358  "42.500835880696 -122.848327688675;42.5008428533324 -122.848245279942" */
1359  a1 = lwgeom_area_sphere(lwg, &s);
1360  CU_ASSERT_DOUBLE_EQUAL(a1, 89.721147136698008, 0.1);
1361  /* spheroid: Planimeter -E -p 20 -r --input-string \
1362  "42.5007249610493 -122.848227067007;42.5007179884263 -122.848309475585;"\
1363  "42.500835880696 -122.848327688675;42.5008428533324 -122.848245279942" */
1364  a2 = lwgeom_area_spheroid(lwg, &s);
1365  CU_ASSERT_DOUBLE_EQUAL(a2, 89.868413479309585, 0.1);
1366  lwgeom_free(lwg);
1367 
1368  /* Big-ass polygon */
1369  lwg = lwgeom_from_wkt("POLYGON((-2 3, -2 4, -1 4, -1 3, -2 3))", LW_PARSER_CHECK_NONE);
1370  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1371  /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE -r --input-string "3 -2;4 -2;4 -1;3 -1" */
1372  a1 = lwgeom_area_sphere(lwg, &s);
1373  CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.106982993974659, 0.1);
1374  /* spheroid: Planimeter -E -p 20 -r --input-string "3 -2;4 -2;4 -1;3 -1" */
1375 #if PROJ_GEODESIC
1376  // printf("XXXXX %d\n", PJ_VERSION);
1377  a2 = lwgeom_area_spheroid(lwg, &s);
1378  CU_ASSERT_DOUBLE_EQUAL(a2, 12286884908.946891319597874, 0.1);
1379 #endif
1380  lwgeom_free(lwg);
1381 
1382  /* One-degree square */
1383  lwg = lwgeom_from_wkt("POLYGON((8.5 2,8.5 1,9.5 1,9.5 2,8.5 2))", LW_PARSER_CHECK_NONE);
1384  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1385  /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE --input-string "2 8.5;1 8.5;1 9.5;2 9.5" */
1386  a1 = lwgeom_area_sphere(lwg, &s);
1387  CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.368023059138681, 0.1);
1388  /* spheroid: Planimeter -E -p 20 --input-string "2 8.5;1 8.5;1 9.5;2 9.5" */
1389 #if PROJ_GEODESIC
1390  a2 = lwgeom_area_spheroid(lwg, &s);
1391  CU_ASSERT_DOUBLE_EQUAL(a2, 12305128751.042900673161556, 0.1);
1392 #endif
1393  lwgeom_free(lwg);
1394 
1395  /* One-degree square *near* the antimeridian */
1396  lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,178.5 1,178.5 2,179.5 2))", LW_PARSER_CHECK_NONE);
1397  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1398  /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE -r --input-string "2 179.5;1 179.5;1 178.5;2 178.5" */
1399  a1 = lwgeom_area_sphere(lwg, &s);
1400  CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.368023059138681, 0.1);
1401  /* spheroid: Planimeter -E -p 20 -r --input-string "2 179.5;1 179.5;1 178.5;2 178.5" */
1402 #if PROJ_GEODESIC
1403  a2 = lwgeom_area_spheroid(lwg, &s);
1404  CU_ASSERT_DOUBLE_EQUAL(a2, 12305128751.042900673161556, 0.1);
1405 #endif
1406  lwgeom_free(lwg);
1407 
1408  /* One-degree square *across* the antimeridian */
1409  lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", LW_PARSER_CHECK_NONE);
1410  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1411  /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE --input-string "2 179.5;1 179.5;1 -179.5;2 -179.5" */
1412  a1 = lwgeom_area_sphere(lwg, &s);
1413  CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.368023059138681, 0.1);
1414  /* spheroid: Planimeter -E -p 20 --input-string "2 179.5;1 179.5;1 -179.5;2 -179.5" */
1415 #if PROJ_GEODESIC
1416  a2 = lwgeom_area_spheroid(lwg, &s);
1417  CU_ASSERT_DOUBLE_EQUAL(a2, 12305128751.042900673161556, 0.1);
1418 #endif
1419  lwgeom_free(lwg);
1420 }
1421 
1422 static void test_gbox_utils(void)
1423 {
1424  LWGEOM *lwg;
1425  GBOX gbox;
1426  double a1, a2;
1427  SPHEROID s;
1428  POINT2D pt;
1429 
1430  /* Init to WGS84 */
1431  spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
1432 
1433  gbox.flags = gflags(0, 0, 1);
1434 
1435  /* One-degree square by equator */
1436  lwg = lwgeom_from_wkt("POLYGON((1 20,1 21,2 21,2 20,1 20))", LW_PARSER_CHECK_NONE);
1437  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1438  a1 = gbox_angular_width(&gbox);
1439  a2 = gbox_angular_height(&gbox);
1440  CU_ASSERT_DOUBLE_EQUAL(a1, 0.0177951, 0.0000001);
1441  CU_ASSERT_DOUBLE_EQUAL(a2, 0.017764, 0.0000001);
1442  lwgeom_free(lwg);
1443 
1444  /* One-degree square *across* antimeridian */
1445  lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", LW_PARSER_CHECK_NONE);
1446  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1447  a1 = gbox_angular_width(&gbox);
1448  a2 = gbox_angular_height(&gbox);
1449  //printf("a1=%g a2=%g\n", a1, a2);
1450  CU_ASSERT_DOUBLE_EQUAL(a1, 0.0174613, 0.0000001);
1451  CU_ASSERT_DOUBLE_EQUAL(a2, 0.0174553, 0.0000001);
1452  lwgeom_free(lwg);
1453 
1454  /* One-degree square *across* antimeridian */
1455  lwg = lwgeom_from_wkt("POLYGON((178.5 2,178.5 1,-179.5 1,-179.5 2,178.5 2))", LW_PARSER_CHECK_NONE);
1456  lwgeom_calculate_gbox_geodetic(lwg, &gbox);
1457  gbox_centroid(&gbox, &pt);
1458  //printf("POINT(%g %g)\n", pt.x, pt.y);
1459  CU_ASSERT_DOUBLE_EQUAL(pt.x, 179.5, 0.0001);
1460  CU_ASSERT_DOUBLE_EQUAL(pt.y, 1.50024, 0.0001);
1461  lwgeom_free(lwg);
1462 
1463 }
1464 
1465 static void test_vector_angle(void)
1466 {
1467  POINT3D p1, p2;
1468  double angle;
1469 
1470  memset(&p1, 0, sizeof(POINT3D));
1471  memset(&p2, 0, sizeof(POINT3D));
1472 
1473  p1.x = 1.0;
1474  p2.y = 1.0;
1475  angle = vector_angle(&p1, &p2);
1476  CU_ASSERT_DOUBLE_EQUAL(angle, M_PI_2, 0.00001);
1477 
1478  p1.x = p2.y = 0.0;
1479  p1.y = 1.0;
1480  p2.x = 1.0;
1481  angle = vector_angle(&p1, &p2);
1482  CU_ASSERT_DOUBLE_EQUAL(angle, M_PI_2, 0.00001);
1483 
1484  p2.y = p2.x = 1.0;
1485  normalize(&p2);
1486  angle = vector_angle(&p1, &p2);
1487  CU_ASSERT_DOUBLE_EQUAL(angle, M_PI_4, 0.00001);
1488 
1489  p2.x = p2.y = p2.z = 1.0;
1490  normalize(&p2);
1491  angle = vector_angle(&p1, &p2);
1492  CU_ASSERT_DOUBLE_EQUAL(angle, 0.955317, 0.00001);
1493  //printf ("angle = %g\n\n", angle);
1494 }
1495 
1496 static void test_vector_rotate(void)
1497 {
1498  POINT3D p1, p2, n;
1499  double angle;
1500 
1501  memset(&p1, 0, sizeof(POINT3D));
1502  memset(&p2, 0, sizeof(POINT3D));
1503  memset(&n, 0, sizeof(POINT3D));
1504 
1505  p1.x = 1.0;
1506  p2.y = 1.0;
1507  angle = M_PI_4;
1508  vector_rotate(&p1, &p2, angle, &n);
1509  //printf("%g %g %g\n\n", n.x, n.y, n.z);
1510  CU_ASSERT_DOUBLE_EQUAL(n.x, 0.707107, 0.00001);
1511 
1512  angle = 2*M_PI/400000000;
1513  vector_rotate(&p1, &p2, angle, &n);
1514  //printf("%.21g %.21g %.21g\n\n", n.x, n.y, n.z);
1515  CU_ASSERT_DOUBLE_EQUAL(n.x, 0.999999999999999888978, 0.0000000000000001);
1516  CU_ASSERT_DOUBLE_EQUAL(n.y, 1.57079632679489654446e-08, 0.0000000000000001);
1517 
1518  angle = 0;
1519  vector_rotate(&p1, &p2, angle, &n);
1520  //printf("%.16g %.16g %.16g\n\n", n.x, n.y, n.z);
1521  CU_ASSERT_DOUBLE_EQUAL(n.x, 1.0, 0.00000001);
1522 }
1523 
1524 static void test_lwgeom_segmentize_sphere(void)
1525 {
1526  LWGEOM *lwg1, *lwg2;
1527  LWLINE *lwl;
1528  double max = 100000.0 / WGS84_RADIUS;
1529  //char *wkt;
1530 
1531  /* Simple case */
1532  lwg1 = lwgeom_from_wkt("LINESTRING(0 20, 5 20)", LW_PARSER_CHECK_NONE);
1533  lwg2 = lwgeom_segmentize_sphere(lwg1, max);
1534  lwl = (LWLINE*)lwg2;
1535  // printf("%s\n", lwgeom_to_ewkt(lwg2));
1536  CU_ASSERT_EQUAL(lwl->points->npoints, 9);
1537  lwgeom_free(lwg1);
1538  lwgeom_free(lwg2);
1539  //lwfree(wkt);
1540 
1541  return;
1542 }
1543 
1544 static void test_lwgeom_area_sphere(void)
1545 {
1546  LWGEOM *lwg;
1547  double area;
1548  SPHEROID s;
1549 
1550  /* Init to WGS84 */
1551  spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
1552 
1553  /* Simple case */
1554  lwg = lwgeom_from_wkt("POLYGON((1 1, 1 2, 2 2, 2 1, 1 1))", LW_PARSER_CHECK_NONE);
1555  area = lwgeom_area_sphere(lwg, &s);
1556 
1557  CU_ASSERT_DOUBLE_EQUAL(area, 12360265021.3561, 1.0);
1558  lwgeom_free(lwg);
1559 
1560  /* Robustness tests, from ticket #3393 */
1561  lwg = lwgeom_from_wkt("POLYGON((0 78.703946026663,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026663))", LW_PARSER_CHECK_NONE);
1562  area = lwgeom_area_sphere(lwg, &s);
1563  CU_ASSERT_DOUBLE_EQUAL(area, 127516467322130, 1.0);
1564  lwgeom_free(lwg);
1565 
1566  lwg = lwgeom_from_wkt("POLYGON((0 78.703946026662,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026662))", LW_PARSER_CHECK_NONE);
1567  area = lwgeom_area_sphere(lwg, &s);
1568  CU_ASSERT_DOUBLE_EQUAL(area, 127516467322130, 1.0);
1569  lwgeom_free(lwg);
1570 
1571  lwg = lwgeom_from_wkt("POLYGON((0 78.703946026664,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026664))", LW_PARSER_CHECK_NONE);
1572  area = lwgeom_area_sphere(lwg, &s);
1573  CU_ASSERT_DOUBLE_EQUAL(area, 127516467322130, 1.0);
1574  lwgeom_free(lwg);
1575  /* end #3393 */
1576 }
1577 
1578 static void test_gbox_to_string_truncated(void)
1579 {
1580  GBOX g = {
1581  .flags = 0,
1582  .xmin = -DBL_MAX,
1583  .xmax = -DBL_MAX,
1584  .ymin = -DBL_MAX,
1585  .ymax = -DBL_MAX,
1586  .zmin = -DBL_MAX,
1587  .zmax = -DBL_MAX,
1588  .mmin = -DBL_MAX,
1589  .mmax = -DBL_MAX,
1590  };
1591  FLAGS_SET_Z(g.flags, 1);
1592  FLAGS_SET_M(g.flags, 1);
1593  char *c = gbox_to_string(&g);
1594 
1595  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))");
1596 
1597  lwfree(c);
1598 }
1599 
1600 /*
1601 ** Used by test harness to register the tests in this file.
1602 */
1603 void geodetic_suite_setup(void);
1604 void geodetic_suite_setup(void)
1605 {
1606  CU_pSuite suite = CU_add_suite("geodetic", NULL, NULL);
1607  PG_ADD_TEST(suite, test_sphere_direction);
1608  PG_ADD_TEST(suite, test_sphere_project);
1609  PG_ADD_TEST(suite, test_lwgeom_area_sphere);
1610  PG_ADD_TEST(suite, test_gbox_from_spherical_coordinates);
1611  PG_ADD_TEST(suite, test_gserialized_get_gbox_geocentric);
1612  PG_ADD_TEST(suite, test_clairaut);
1613  PG_ADD_TEST(suite, test_edge_intersection);
1614  PG_ADD_TEST(suite, test_edge_intersects);
1615  PG_ADD_TEST(suite, test_edge_distance_to_point);
1616  PG_ADD_TEST(suite, test_edge_distance_to_edge);
1617  PG_ADD_TEST(suite, test_lwgeom_distance_sphere);
1618  PG_ADD_TEST(suite, test_lwgeom_check_geodetic);
1619  PG_ADD_TEST(suite, test_gserialized_from_lwgeom);
1620  PG_ADD_TEST(suite, test_spheroid_distance);
1621  PG_ADD_TEST(suite, test_spheroid_area);
1622  PG_ADD_TEST(suite, test_lwpoly_covers_point2d);
1623  PG_ADD_TEST(suite, test_gbox_utils);
1624  PG_ADD_TEST(suite, test_vector_angle);
1625  PG_ADD_TEST(suite, test_vector_rotate);
1626  PG_ADD_TEST(suite, test_lwgeom_segmentize_sphere);
1627  PG_ADD_TEST(suite, test_ptarray_contains_point_sphere);
1628  PG_ADD_TEST(suite, test_ptarray_contains_point_sphere_iowa);
1629  PG_ADD_TEST(suite, test_gbox_to_string_truncated);
1630 }
test_gserialized_from_lwgeom
static void test_gserialized_from_lwgeom(void)
Definition: cu_geodetic.c:894
test_ptarray_contains_point_sphere
static void test_ptarray_contains_point_sphere(void)
Definition: cu_geodetic.c:933
test_spheroid_distance
static void test_spheroid_distance(void)
Definition: cu_geodetic.c:1275
point_rad2deg
static void point_rad2deg(GEOGRAPHIC_POINT *p)
Convert a point from radians to degrees.
Definition: cu_geodetic.c:63
test_clairaut
static void test_clairaut(void)
Definition: cu_geodetic.c:363
test_sphere_direction
static void test_sphere_direction(void)
Definition: cu_geodetic.c:69
edge_deg2rad
static void edge_deg2rad(GEOGRAPHIC_EDGE *e)
Convert an edge from degrees to radians.
Definition: cu_geodetic.c:32
test_lwgeom_segmentize_sphere
static void test_lwgeom_segmentize_sphere(void)
Definition: cu_geodetic.c:1524
test_gbox_from_spherical_coordinates
static void test_gbox_from_spherical_coordinates(void)
Definition: cu_geodetic.c:224
test_vector_rotate
static void test_vector_rotate(void)
Definition: cu_geodetic.c:1496
point_set
static void point_set(double lon, double lat, GEOGRAPHIC_POINT *p)
Definition: cu_geodetic.c:356
test_gbox_to_string_truncated
static void test_gbox_to_string_truncated(void)
Definition: cu_geodetic.c:1578
test_lwpoly_covers_point2d
static void test_lwpoly_covers_point2d(void)
Definition: cu_geodetic.c:1122
test_sphere_project
static void test_sphere_project(void)
Definition: cu_geodetic.c:92
test_edge_distance_to_edge
static void test_edge_distance_to_edge(void)
Definition: cu_geodetic.c:767
edge_set
static void edge_set(double lon1, double lat1, double lon2, double lat2, GEOGRAPHIC_EDGE *e)
Definition: cu_geodetic.c:347
point_deg2rad
static void point_deg2rad(GEOGRAPHIC_POINT *p)
Convert an edge from radians to degrees.
Definition: cu_geodetic.c:54
test_edge_intersects
static void test_edge_intersects(void)
Definition: cu_geodetic.c:583
geodetic_suite_setup
void geodetic_suite_setup(void)
Definition: cu_geodetic.c:1604
lwgeom_over_gserialized
static LWGEOM * lwgeom_over_gserialized(char *wkt, GSERIALIZED **g)
Definition: cu_geodetic.c:795
test_edge_distance_to_point
static void test_edge_distance_to_point(void)
Definition: cu_geodetic.c:730
test_lwgeom_area_sphere
static void test_lwgeom_area_sphere(void)
Definition: cu_geodetic.c:1544
RANDOM_TEST
#define RANDOM_TEST
Definition: cu_geodetic.c:27
test_ptarray_contains_point_sphere_iowa
static void test_ptarray_contains_point_sphere_iowa(void)
Definition: cu_geodetic.c:1157
test_edge_intersection
static void test_edge_intersection(void)
Definition: cu_geodetic.c:401
test_lwgeom_check_geodetic
static void test_lwgeom_check_geodetic(void)
Definition: cu_geodetic.c:806
line2pts
static void line2pts(const char *wkt, POINT3D *A1, POINT3D *A2)
Definition: cu_geodetic.c:561
test_spheroid_area
static void test_spheroid_area(void)
Definition: cu_geodetic.c:1341
test_gbox_utils
static void test_gbox_utils(void)
Definition: cu_geodetic.c:1422
test_gserialized_get_gbox_geocentric
static void test_gserialized_get_gbox_geocentric(void)
Definition: cu_geodetic.c:295
test_vector_angle
static void test_vector_angle(void)
Definition: cu_geodetic.c:1465
test_lwgeom_distance_sphere
static void test_lwgeom_distance_sphere(void)
Definition: cu_geodetic.c:1177
gbox_data_length
int gbox_data_length
Definition: cu_geodetic_data.h:12
iowa_data
char * iowa_data
Definition: cu_geodetic_data.h:74
gbox_data
char gbox_data[][512]
Definition: cu_geodetic_data.h:13
cu_geodetic_data.h
s
char * s
Definition: cu_in_wkt.c:23
gbox_to_string
char * gbox_to_string(const GBOX *gbox)
Allocate a string representation of the GBOX, based on dimensionality of flags.
Definition: g_box.c:399
lwgeom_from_gserialized
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Definition: g_serialized.c:1530
gserialized_get_type
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
gserialized_from_lwgeom
GSERIALIZED * gserialized_from_lwgeom(LWGEOM *geom, size_t *size)
Allocate a new GSERIALIZED from an LWGEOM.
Definition: g_serialized.c:1176
gflags
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:145
PG_ADD_TEST
#define PG_ADD_TEST(suite, testfunc)
Definition: liblwgeom/cunit/cu_tester.h:20
ASSERT_STRING_EQUAL
#define ASSERT_STRING_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:43
ptarray_construct_reference_data
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
WGS84_MAJOR_AXIS
#define WGS84_MAJOR_AXIS
Definition: liblwgeom.h:128
lwline_as_lwgeom
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:330
LW_FALSE
#define LW_FALSE
Definition: liblwgeom.h:77
lwgeom_segmentize_sphere
LWGEOM * lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length)
Derive a new geometry with vertices added to ensure no vertex is more than max_seg_length (in radians...
Definition: lwgeodetic.c:1743
lwgeom_from_hexwkb
LWGEOM * lwgeom_from_hexwkb(const char *hexwkb, const char check)
Definition: lwin_wkb.c:809
lwgeom_distance_spheroid
double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance)
Calculate the geodetic distance from lwgeom1 to lwgeom2 on the spheroid.
Definition: lwgeodetic.c:2187
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
LW_PARSER_CHECK_NONE
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:2005
WGS84_MINOR_AXIS
#define WGS84_MINOR_AXIS
Definition: liblwgeom.h:130
MULTILINETYPE
#define MULTILINETYPE
Definition: liblwgeom.h:89
lwgeom_calculate_gbox_geodetic
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:3028
WGS84_RADIUS
#define WGS84_RADIUS
Definition: liblwgeom.h:131
lwgeom_area_spheroid
double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the geodetic area of a lwgeom on the spheroid.
Definition: lwspheroid.c:647
lwgeom_check_geodetic
int lwgeom_check_geodetic(const LWGEOM *geom)
Check that coordinates of LWGEOM are all within the geodetic range (-180, -90, 180,...
Definition: lwgeodetic.c:3131
getPoint2d_p
int getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
Definition: lwgeom_api.c:348
POINTTYPE
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:85
lwfree
void lwfree(void *mem)
Definition: lwutil.c:244
POLYGONTYPE
#define POLYGONTYPE
Definition: liblwgeom.h:87
spheroid_init
void spheroid_init(SPHEROID *s, double a, double b)
Initialize a spheroid object for use in geodetic functions.
Definition: lwspheroid.c:39
lwgeom_calculate_gbox
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
lwgeom_from_wkt
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition: lwin_wkt.c:904
FLAGS_SET_GEODETIC
#define FLAGS_SET_GEODETIC(flags, value)
Definition: liblwgeom.h:149
LW_TRUE
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:76
FLAGS_SET_M
#define FLAGS_SET_M(flags, value)
Definition: liblwgeom.h:147
SRID_UNKNOWN
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:188
FLAGS_SET_Z
#define FLAGS_SET_Z(flags, value)
Definition: liblwgeom.h:146
lwline_construct
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
lwline_free
void lwline_free(LWLINE *line)
Definition: lwline.c:76
lwgeom_area_sphere
double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the geodetic area of a lwgeom on the sphere.
Definition: lwgeodetic.c:2031
gbox_angular_height
double gbox_angular_height(const GBOX *gbox)
GBOX utility functions to figure out coverage/location on the globe.
Definition: lwgeodetic.c:188
gbox_angular_width
double gbox_angular_width(const GBOX *gbox)
Returns the angular width (longitudinal span) of the box in radians.
Definition: lwgeodetic.c:215
gbox_centroid
int gbox_centroid(const GBOX *gbox, POINT2D *out)
Computes the average(ish) center of the box and returns success.
Definition: lwgeodetic.c:267
liblwgeom_internal.h
clairaut_geographic
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:1101
lwpoly_covers_point2d
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:2518
normalize
void normalize(POINT3D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:615
longitude_radians_normalize
double longitude_radians_normalize(double lon)
Convert a longitude to the range of -PI,PI.
Definition: lwgeodetic.c:50
clairaut_cartesian
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:1076
vector_rotate
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:573
ptarray_contains_point_sphere
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:3645
latitude_radians_normalize
double latitude_radians_normalize(double lat)
Convert a latitude to the range of -PI/2,PI/2.
Definition: lwgeodetic.c:78
sphere_project
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:1314
ll2cart
void ll2cart(const POINT2D *g, POINT3D *p)
Convert lon/lat coordinates to cartesian coordinates on unit sphere.
Definition: lwgeodetic.c:423
gbox_geocentric_slow
int gbox_geocentric_slow
For testing geodetic bounding box, we have a magic global variable.
Definition: lwgeodetic.c:36
robust_cross_product
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:634
geographic_point_init
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:180
unit_normal
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:541
sphere_distance
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:948
sphere_direction
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:975
edge_intersection
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:1127
edge_intersects
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:3540
edge_distance_to_point
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
Definition: lwgeodetic.c:1218
geog2cart
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesian coordinates on unit sphere.
Definition: lwgeodetic.c:404
edge_distance_to_edge
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:1269
vector_angle
double vector_angle(const POINT3D *v1, const POINT3D *v2)
Angle between two unit vectors.
Definition: lwgeodetic.c:505
rad2deg
#define rad2deg(r)
Definition: lwgeodetic.h:80
spheroid_distance
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:191
PIR_COLINEAR
#define PIR_COLINEAR
Definition: lwgeodetic.h:88
PIR_INTERSECTS
#define PIR_INTERSECTS
Definition: lwgeodetic.h:87
deg2rad
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:79
PIR_A_TOUCH_RIGHT
#define PIR_A_TOUCH_RIGHT
Definition: lwgeodetic.h:89
PIR_B_TOUCH_RIGHT
#define PIR_B_TOUCH_RIGHT
Definition: lwgeodetic.h:91
PIR_B_TOUCH_LEFT
#define PIR_B_TOUCH_LEFT
Definition: lwgeodetic.h:92
area
Datum area(PG_FUNCTION_ARGS)
Definition: lwgeom_backend_api.c:193
ovdump.type
type
Definition: ovdump.py:41
cu_tester.h
GBOX::ymax
double ymax
Definition: liblwgeom.h:298
GBOX::zmax
double zmax
Definition: liblwgeom.h:300
GBOX::xmax
double xmax
Definition: liblwgeom.h:296
GBOX::zmin
double zmin
Definition: liblwgeom.h:299
GBOX::ymin
double ymin
Definition: liblwgeom.h:297
GBOX::xmin
double xmin
Definition: liblwgeom.h:295
GBOX::flags
uint8_t flags
Definition: liblwgeom.h:294
GEOGRAPHIC_EDGE::start
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:63
GEOGRAPHIC_EDGE::end
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:64
GEOGRAPHIC_EDGE
Two-point great circle segment from a to b.
Definition: lwgeodetic.h:62
GEOGRAPHIC_POINT::lat
double lat
Definition: lwgeodetic.h:55
GEOGRAPHIC_POINT::lon
double lon
Definition: lwgeodetic.h:54
GEOGRAPHIC_POINT
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:53
LWGEOM::flags
uint8_t flags
Definition: liblwgeom.h:400
LWLINE::points
POINTARRAY * points
Definition: liblwgeom.h:425
LWPOLY::rings
POINTARRAY ** rings
Definition: liblwgeom.h:460
POINT2D::y
double y
Definition: liblwgeom.h:331
POINT2D::x
double x
Definition: liblwgeom.h:331
POINT3D::z
double z
Definition: liblwgeom.h:343
POINT3D::x
double x
Definition: liblwgeom.h:343
POINT3D::y
double y
Definition: liblwgeom.h:343
POINTARRAY::npoints
uint32_t npoints
Definition: liblwgeom.h:374
uint32_t
unsigned int uint32_t
Definition: uthash.h:78
uint8_t
unsigned char uint8_t
Definition: uthash.h:79