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
29 /**
30 * Convert an edge from degrees to radians.
31 */
edge_deg2rad(GEOGRAPHIC_EDGE * e)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
40 /**
41 * Convert an edge from radians to degrees.
42 static void edge_rad2deg(GEOGRAPHIC_EDGE *e)
43 {
44 (e->start).lat = rad2deg((e->start).lat);
45 (e->end).lat = rad2deg((e->end).lat);
46 (e->start).lon = rad2deg((e->start).lon);
47 (e->end).lon = rad2deg((e->end).lon);
48 }
49 */
50
51 /**
52 * Convert a point from degrees to radians.
53 */
point_deg2rad(GEOGRAPHIC_POINT * p)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
60 /**
61 * Convert a point from radians to degrees.
62 */
point_rad2deg(GEOGRAPHIC_POINT * p)63 static void point_rad2deg(GEOGRAPHIC_POINT *p)
64 {
65 p->lat = rad2deg(p->lat);
66 p->lon = rad2deg(p->lon);
67 }
68
test_sphere_direction(void)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
test_sphere_project(void)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
173 /**
174 * Tests the relative numerical stability of the "robust" and
175 * naive cross product calculation methods.
176 */
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
test_gbox_from_spherical_coordinates(void)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
test_gserialized_get_gbox_geocentric(void)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
edge_set(double lon1,double lat1,double lon2,double lat2,GEOGRAPHIC_EDGE * e)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
point_set(double lon,double lat,GEOGRAPHIC_POINT * p)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
test_clairaut(void)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
test_edge_intersection(void)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
line2pts(const char * wkt,POINT3D * A1,POINT3D * A2)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
test_edge_intersects(void)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
test_edge_distance_to_point(void)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
test_edge_distance_to_edge(void)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 */
lwgeom_over_gserialized(char * wkt,GSERIALIZED ** g)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
test_lwgeom_check_geodetic(void)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
test_gserialized_from_lwgeom(void)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
test_ptarray_contains_point_sphere(void)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
test_lwpoly_covers_point2d(void)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
test_ptarray_contains_point_sphere_iowa(void)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
test_lwgeom_distance_sphere(void)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
test_spheroid_distance(void)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
test_spheroid_area(void)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
test_gbox_utils(void)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
test_vector_angle(void)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
test_vector_rotate(void)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
test_lwgeom_segmentize_sphere(void)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
test_lwgeom_area_sphere(void)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
test_gbox_to_string_truncated(void)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);
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 }
1631