source: trunk/liblwgeom/cunit/cu_geodetic.c @ 5345

Last change on this file since 5345 was 5345, checked in by pramsey, 7 years ago

Fix potential corner case in sphere area calculation (#451)

  • Property svn:keywords set to Author Date Id Revision
File size: 34.2 KB
Line 
1/**********************************************************************
2 * $Id: cu_geodetic.c 5345 2010-02-25 15:13:48Z pramsey $
3 *
4 * PostGIS - Spatial Types for PostgreSQL
5 * http://postgis.refractions.net
6 * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
7 *
8 * This is free software; you can redistribute and/or modify it under
9 * the terms of the GNU General Public Licence. See the COPYING file.
10 *
11 **********************************************************************/
12
13#include <stdio.h>
14#include <stdlib.h>
15#include <string.h>
16#include "CUnit/Basic.h"
17
18#include "lwgeodetic.h"
19#include "cu_tester.h"
20
21#define RANDOM_TEST 0
22
23/**
24* Convert an edge from degrees to radians.
25*/
26static void edge_deg2rad(GEOGRAPHIC_EDGE *e)
27{
28        (e->start).lat = deg2rad((e->start).lat);
29        (e->end).lat = deg2rad((e->end).lat);
30        (e->start).lon = deg2rad((e->start).lon);
31        (e->end).lon = deg2rad((e->end).lon);
32}
33
34/**
35* Convert an edge from radians to degrees.
36static void edge_rad2deg(GEOGRAPHIC_EDGE *e)
37{
38        (e->start).lat = rad2deg((e->start).lat);
39        (e->end).lat = rad2deg((e->end).lat);
40        (e->start).lon = rad2deg((e->start).lon);
41        (e->end).lon = rad2deg((e->end).lon);
42}
43*/
44
45/**
46* Convert a point from degrees to radians.
47*/
48static void point_deg2rad(GEOGRAPHIC_POINT *p)
49{
50        p->lat = latitude_radians_normalize(deg2rad(p->lat));
51        p->lon = longitude_radians_normalize(deg2rad(p->lon));
52}
53
54/**
55* Convert a point from radians to degrees.
56*/
57static void point_rad2deg(GEOGRAPHIC_POINT *p)
58{
59        p->lat = rad2deg(p->lat);
60        p->lon = rad2deg(p->lon);
61}
62
63static void test_signum(void)
64{
65        CU_ASSERT_EQUAL(signum(-5.0),-1);
66        CU_ASSERT_EQUAL(signum(5.0),1);
67}
68
69
70static void test_gbox_from_spherical_coordinates(void)
71{
72#if RANDOM_TEST
73        const double gtolerance = 0.000001;
74        const int loops = RANDOM_TEST;
75        int i;
76        double ll[64];
77        GBOX *gbox;
78        GBOX *gbox_slow;
79        int rndlat;
80        int rndlon;
81
82        POINTARRAY *pa;
83        LWLINE *lwline;
84        GSERIALIZED *g;
85
86        ll[0] = -3.083333333333333333333333333333333;
87        ll[1] = 9.83333333333333333333333333333333;
88        ll[2] = 15.5;
89        ll[3] = -5.25;
90
91        pa = pointArray_construct((uchar*)ll, 0, 0, 2);
92        lwline = lwline_construct(-1, 0, pa);
93
94        srandomdev();
95
96        for ( i = 0; i < loops; i++ )
97        {
98                rndlat = (int)(90.0 - 180.0 * (double)random() / pow(2.0, 31.0));
99                rndlon = (int)(180.0 - 360.0 * (double)random() / pow(2.0, 31.0));
100                ll[0] = (double)rndlon;
101                ll[1] = (double)rndlat;
102
103                rndlat = (int)(90.0 - 180.0 * (double)random() / pow(2.0, 31.0));
104                rndlon = (int)(180.0 - 360.0 * (double)random() / pow(2.0, 31.0));
105                ll[2] = (double)rndlon;
106                ll[3] = (double)rndlat;
107
108                g = gserialized_from_lwgeom((LWGEOM*)lwline, 1, 0);
109                FLAGS_SET_GEODETIC(g->flags, 1);
110                gbox_geocentric_slow = LW_FALSE;
111                gbox = gserialized_calculate_gbox_geocentric(g);
112                gbox_geocentric_slow = LW_TRUE;
113                gbox_slow = gserialized_calculate_gbox_geocentric(g);
114                gbox_geocentric_slow = LW_FALSE;
115                lwfree(g);
116
117                if (
118                    ( fabs( gbox->xmin - gbox_slow->xmin ) > gtolerance ) ||
119                    ( fabs( gbox->xmax - gbox_slow->xmax ) > gtolerance ) ||
120                    ( fabs( gbox->ymin - gbox_slow->ymin ) > gtolerance ) ||
121                    ( fabs( gbox->ymax - gbox_slow->ymax ) > gtolerance ) ||
122                    ( fabs( gbox->zmin - gbox_slow->zmin ) > gtolerance ) ||
123                    ( fabs( gbox->zmax - gbox_slow->zmax ) > gtolerance ) )
124                {
125                        printf("\n-------\n");
126                        printf("If you are seeing this, cut and paste it, it is a randomly generated test case!\n");
127                        printf("LOOP: %d\n", i);
128                        printf("SEGMENT (Lon Lat): (%.9g %.9g) (%.9g %.9g)\n", ll[0], ll[1], ll[2], ll[3]);
129                        printf("CALC: %s\n", gbox_to_string(gbox));
130                        printf("SLOW: %s\n", gbox_to_string(gbox_slow));
131                        printf("-------\n\n");
132                        CU_FAIL_FATAL(Slow (GOOD) and fast (CALC) box calculations returned different values!!);
133                }
134
135                lwfree(gbox);
136                lwfree(gbox_slow);
137        }
138
139        lwfree(lwline);
140        lwfree(pa);
141#endif /* RANDOM_TEST */
142}
143
144#include "cu_geodetic_data.h"
145
146static void test_gserialized_get_gbox_geocentric(void)
147{
148        LWGEOM *lwg;
149        GSERIALIZED *g;
150        GBOX *gbox, *gbox_slow;
151        int i;
152
153        for ( i = 0; i < gbox_data_length; i++ )
154        {
155#if 0
156//              if ( i != 0 ) continue; /* skip our bad case */
157                printf("\n\n------------\n");
158                printf("%s\n", gbox_data[i]);
159#endif
160                lwg = lwgeom_from_ewkt(gbox_data[i], PARSER_CHECK_NONE);
161                g = gserialized_from_lwgeom(lwg, 1, 0);
162                FLAGS_SET_GEODETIC(g->flags, 1);
163                lwgeom_free(lwg);
164                gbox_geocentric_slow = LW_FALSE;
165                gbox = gserialized_calculate_gbox_geocentric(g);
166                gbox_geocentric_slow = LW_TRUE;
167                gbox_slow = gserialized_calculate_gbox_geocentric(g);
168                gbox_geocentric_slow = LW_FALSE;
169#if 0
170                printf("\nCALC: %s\n", gbox_to_string(gbox));
171                printf("GOOD: %s\n", gbox_to_string(gbox_slow));
172                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));
173                printf("------------\n");
174#endif
175                CU_ASSERT_DOUBLE_EQUAL(gbox->xmin, gbox_slow->xmin, 0.000001);
176                CU_ASSERT_DOUBLE_EQUAL(gbox->ymin, gbox_slow->ymin, 0.000001);
177                CU_ASSERT_DOUBLE_EQUAL(gbox->zmin, gbox_slow->zmin, 0.000001);
178                CU_ASSERT_DOUBLE_EQUAL(gbox->xmax, gbox_slow->xmax, 0.000001);
179                CU_ASSERT_DOUBLE_EQUAL(gbox->ymax, gbox_slow->ymax, 0.000001);
180                CU_ASSERT_DOUBLE_EQUAL(gbox->zmax, gbox_slow->zmax, 0.000001);
181                lwfree(g);
182                lwfree(gbox);
183                lwfree(gbox_slow);
184        }
185
186}
187
188/*
189* Build LWGEOM on top of *aligned* structure so we can use the read-only
190* point access methods on them.
191static LWGEOM* lwgeom_over_gserialized(char *wkt)
192{
193        LWGEOM *lwg;
194        GSERIALIZED *g;
195
196        lwg = lwgeom_from_ewkt(wkt, PARSER_CHECK_NONE);
197        g = gserialized_from_lwgeom(lwg, 1, 0);
198        lwgeom_free(lwg);
199        return lwgeom_from_gserialized(g);
200}
201*/
202
203static void edge_set(double lon1, double lat1, double lon2, double lat2, GEOGRAPHIC_EDGE *e)
204{
205        e->start.lon = lon1;
206        e->start.lat = lat1;
207        e->end.lon = lon2;
208        e->end.lat = lat2;
209        edge_deg2rad(e);
210}
211
212static void point_set(double lon, double lat, GEOGRAPHIC_POINT *p)
213{
214        p->lon = lon;
215        p->lat = lat;
216        point_deg2rad(p);
217}
218
219static void test_clairaut(void)
220{
221
222        GEOGRAPHIC_POINT gs, ge;
223        POINT3D vs, ve;
224        GEOGRAPHIC_POINT g_out_top, g_out_bottom, v_out_top, v_out_bottom;
225
226        point_set(-45.0, 60.0, &gs);
227        point_set(135.0, 60.0, &ge);
228
229        geog2cart(&gs, &vs);
230        geog2cart(&ge, &ve);
231
232        clairaut_cartesian(&vs, &ve, &v_out_top, &v_out_bottom);
233        clairaut_geographic(&gs, &ge, &g_out_top, &g_out_bottom);
234
235        CU_ASSERT_DOUBLE_EQUAL(v_out_top.lat, g_out_top.lat, 0.000001);
236        CU_ASSERT_DOUBLE_EQUAL(v_out_top.lon, g_out_top.lon, 0.000001);
237        CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lat, g_out_bottom.lat, 0.000001);
238        CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lon, g_out_bottom.lon, 0.000001);
239
240        gs.lat = 1.3021240033804449;
241        ge.lat = 1.3021240033804449;
242        gs.lon = -1.3387392931438733;
243        ge.lon = 1.80285336044592;
244
245        geog2cart(&gs, &vs);
246        geog2cart(&ge, &ve);
247
248        clairaut_cartesian(&vs, &ve, &v_out_top, &v_out_bottom);
249        clairaut_geographic(&gs, &ge, &g_out_top, &g_out_bottom);
250
251        CU_ASSERT_DOUBLE_EQUAL(v_out_top.lat, g_out_top.lat, 0.000001);
252        CU_ASSERT_DOUBLE_EQUAL(v_out_top.lon, g_out_top.lon, 0.000001);
253        CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lat, g_out_bottom.lat, 0.000001);
254        CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lon, g_out_bottom.lon, 0.000001);
255}
256
257static void test_edge_intersection(void)
258{
259        GEOGRAPHIC_EDGE e1, e2;
260        GEOGRAPHIC_POINT g;
261        int rv;
262
263        /* Covers case, end-to-end intersection */
264        edge_set(50, -10.999999999999998224, -10.0, 50.0, &e1);
265        edge_set(-10.0, 50.0, -10.272779983831613393, -16.937003313332997578, &e2);
266        rv = edge_intersection(&e1, &e2, &g);
267        CU_ASSERT_EQUAL(rv, LW_TRUE);
268
269        /* Medford case, very short segment vs very long one */
270        e1.start.lat = 0.74123572595649878103;
271        e1.start.lon = -2.1496353191142714145;
272        e1.end.lat = 0.74123631950116664058;
273        e1.end.lon = -2.1496353248304860273;
274        e2.start.lat = 0.73856343764436815924;
275        e2.start.lon = -2.1461493501950630325;
276        e2.end.lat = 0.70971354024834598651;
277        e2.end.lon = 2.1082194552519770703;
278        rv = edge_intersection(&e1, &e2, &g);
279        CU_ASSERT_EQUAL(rv, LW_FALSE);
280
281        /* Again, this time with a less exact input edge. */
282        edge_set(-123.165031277506, 42.4696787216231, -123.165031605021, 42.4697127292275, &e1);
283        rv = edge_intersection(&e1, &e2, &g);
284        CU_ASSERT_EQUAL(rv, LW_FALSE);
285
286        /* Second Medford case, very short segment vs very long one
287        e1.start.lat = 0.73826546728290887156;
288        e1.start.lon = -2.14426380171833042;
289        e1.end.lat = 0.73826545883786642843;
290        e1.end.lon = -2.1442638997530165668;
291        e2.start.lat = 0.73775469118192538165;
292        e2.start.lon = -2.1436035534281718817;
293        e2.end.lat = 0.71021099548296817705;
294        e2.end.lon = 2.1065275171200439353;
295        rv = edge_intersection(e1, e2, &g);
296        CU_ASSERT_EQUAL(rv, LW_FALSE);
297        */
298
299        /* Intersection at (0 0) */
300        edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
301        edge_set(0.0, -1.0, 0.0, 1.0, &e2);
302        rv = edge_intersection(&e1, &e2, &g);
303        point_rad2deg(&g);
304        CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
305        CU_ASSERT_DOUBLE_EQUAL(g.lon, 0.0, 0.00001);
306        CU_ASSERT_EQUAL(rv, LW_TRUE);
307
308        /*  No intersection at (0 0) */
309        edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
310        edge_set(0.0, -1.0, 0.0, -2.0, &e2);
311        rv = edge_intersection(&e1, &e2, &g);
312        CU_ASSERT_EQUAL(rv, LW_FALSE);
313
314        /*  End touches middle of segment at (0 0) */
315        edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
316        edge_set(0.0, -1.0, 0.0, 0.0, &e2);
317        rv = edge_intersection(&e1, &e2, &g);
318        point_rad2deg(&g);
319#if 0
320        printf("\n");
321        printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon,  e1.start.lat, e1.end.lon,  e1.end.lat);
322        printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon,  e2.start.lat, e2.end.lon,  e2.end.lat);
323        printf("g = (%.15g %.15g)\n", g.lon, g.lat);
324        printf("rv = %d\n", rv);
325#endif
326        CU_ASSERT_DOUBLE_EQUAL(g.lon, 0.0, 0.00001);
327        CU_ASSERT_EQUAL(rv, LW_TRUE);
328
329        /*  End touches end of segment at (0 0) */
330        edge_set(0.0, 0.0, 1.0, 0.0, &e1);
331        edge_set(0.0, -1.0, 0.0, 0.0, &e2);
332        rv = edge_intersection(&e1, &e2, &g);
333        point_rad2deg(&g);
334#if 0
335        printf("\n");
336        printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon,  e1.start.lat, e1.end.lon,  e1.end.lat);
337        printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon,  e2.start.lat, e2.end.lon,  e2.end.lat);
338        printf("g = (%.15g %.15g)\n", g.lon, g.lat);
339        printf("rv = %d\n", rv);
340#endif
341        CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
342        CU_ASSERT_DOUBLE_EQUAL(g.lon, 0.0, 0.00001);
343        CU_ASSERT_EQUAL(rv, LW_TRUE);
344
345        /* Intersection at (180 0) */
346        edge_set(-179.0, -1.0, 179.0, 1.0, &e1);
347        edge_set(-179.0, 1.0, 179.0, -1.0, &e2);
348        rv = edge_intersection(&e1, &e2, &g);
349        point_rad2deg(&g);
350        CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
351        CU_ASSERT_DOUBLE_EQUAL(fabs(g.lon), 180.0, 0.00001);
352        CU_ASSERT_EQUAL(rv, LW_TRUE);
353
354        /* Intersection at (180 0) */
355        edge_set(-170.0, 0.0, 170.0, 0.0, &e1);
356        edge_set(180.0, -10.0, 180.0, 10.0, &e2);
357        rv = edge_intersection(&e1, &e2, &g);
358        point_rad2deg(&g);
359        CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
360        CU_ASSERT_DOUBLE_EQUAL(fabs(g.lon), 180.0, 0.00001);
361        CU_ASSERT_EQUAL(rv, LW_TRUE);
362
363        /* Intersection at north pole */
364        edge_set(-180.0, 80.0, 0.0, 80.0, &e1);
365        edge_set(90.0, 80.0, -90.0, 80.0, &e2);
366        rv = edge_intersection(&e1, &e2, &g);
367        point_rad2deg(&g);
368        CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
369        CU_ASSERT_EQUAL(rv, LW_TRUE);
370
371        /* Equal edges return true */
372        edge_set(45.0, 10.0, 50.0, 20.0, &e1);
373        edge_set(45.0, 10.0, 50.0, 20.0, &e2);
374        rv = edge_intersection(&e1, &e2, &g);
375        point_rad2deg(&g);
376        CU_ASSERT_EQUAL(rv, LW_TRUE);
377
378        /* Parallel edges (same great circle, different end points) return true  */
379        edge_set(40.0, 0.0, 70.0, 0.0, &e1);
380        edge_set(60.0, 0.0, 50.0, 0.0, &e2);
381        rv = edge_intersection(&e1, &e2, &g);
382        point_rad2deg(&g);
383        CU_ASSERT_EQUAL(rv, 2); /* Hack, returning 2 as the 'co-linear' value */
384
385        /* End touches arc at north pole */
386        edge_set(-180.0, 80.0, 0.0, 80.0, &e1);
387        edge_set(90.0, 80.0, -90.0, 90.0, &e2);
388        rv = edge_intersection(&e1, &e2, &g);
389        point_rad2deg(&g);
390#if 0
391        printf("\n");
392        printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon,  e1.start.lat, e1.end.lon,  e1.end.lat);
393        printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon,  e2.start.lat, e2.end.lon,  e2.end.lat);
394        printf("g = (%.15g %.15g)\n", g.lon, g.lat);
395        printf("rv = %d\n", rv);
396#endif
397        CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
398        CU_ASSERT_EQUAL(rv, LW_TRUE);
399
400        /* End touches end at north pole */
401        edge_set(-180.0, 80.0, 0.0, 90.0, &e1);
402        edge_set(90.0, 80.0, -90.0, 90.0, &e2);
403        rv = edge_intersection(&e1, &e2, &g);
404        point_rad2deg(&g);
405#if 0
406        printf("\n");
407        printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon,  e1.start.lat, e1.end.lon,  e1.end.lat);
408        printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon,  e2.start.lat, e2.end.lon,  e2.end.lat);
409        printf("g = (%.15g %.15g)\n", g.lon, g.lat);
410        printf("rv = %d\n", rv);
411#endif
412        CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
413        CU_ASSERT_EQUAL(rv, LW_TRUE);
414
415}
416
417static void test_edge_distance_to_point(void)
418{
419        GEOGRAPHIC_EDGE e;
420        GEOGRAPHIC_POINT g;
421        GEOGRAPHIC_POINT closest;
422        double d;
423
424        /* closest point at origin, one degree away */
425        edge_set(-50.0, 0.0, 50.0, 0.0, &e);
426        point_set(0.0, 1.0, &g);
427        d = edge_distance_to_point(&e, &g, 0);
428        CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001);
429
430        /* closest point at origin, one degree away */
431        edge_set(-50.0, 0.0, 50.0, 0.0, &e);
432        point_set(0.0, 2.0, &g);
433        d = edge_distance_to_point(&e, &g, &closest);
434#if 0
435        printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e.start.lon,  e.start.lat, e.end.lon,  e.end.lat);
436        printf("POINT(%.9g %.9g)\n", g.lon, g.lat);
437        printf("\nDISTANCE == %.8g\n", d);
438#endif
439        CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 90.0, 0.00001);
440        CU_ASSERT_DOUBLE_EQUAL(closest.lat, 0.0, 0.00001);
441        CU_ASSERT_DOUBLE_EQUAL(closest.lon, 0.0, 0.00001);
442
443}
444
445static void test_edge_distance_to_edge(void)
446{
447        GEOGRAPHIC_EDGE e1, e2;
448        GEOGRAPHIC_POINT c1, c2;
449        double d;
450
451        /* closest point at origin, one degree away */
452        edge_set(-50.0, 0.0, 50.0, 0.0, &e1);
453        edge_set(-5.0, 20.0, 0.0, 1.0, &e2);
454        d = edge_distance_to_edge(&e1, &e2, &c1, &c2);
455#if 0
456        printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e1.start.lon,  e1.start.lat, e1.end.lon,  e1.end.lat);
457        printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e2.start.lon,  e2.start.lat, e2.end.lon,  e2.end.lat);
458        printf("\nDISTANCE == %.8g\n", d);
459#endif
460        CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001);
461        CU_ASSERT_DOUBLE_EQUAL(c1.lat, 0.0, 0.00001);
462        CU_ASSERT_DOUBLE_EQUAL(c2.lat, M_PI / 180.0, 0.00001);
463        CU_ASSERT_DOUBLE_EQUAL(c1.lon, 0.0, 0.00001);
464        CU_ASSERT_DOUBLE_EQUAL(c2.lon, 0.0, 0.00001);
465}
466
467
468
469/*
470* Build LWGEOM on top of *aligned* structure so we can use the read-only
471* point access methods on them.
472*/
473static LWGEOM* lwgeom_over_gserialized(char *wkt)
474{
475        LWGEOM *lwg;
476        GSERIALIZED *g;
477
478        lwg = lwgeom_from_ewkt(wkt, PARSER_CHECK_NONE);
479        g = gserialized_from_lwgeom(lwg, 1, 0);
480        lwgeom_free(lwg);
481        return lwgeom_from_gserialized(g);
482}
483
484static void test_lwgeom_check_geodetic(void)
485{
486        LWGEOM *geom;
487        int i = 0;
488
489        char ewkt[][512] =
490        {
491                "POINT(0 0.2)",
492                "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
493                "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))",
494                "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
495                "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)))",
496                "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))))",
497                "POINT(0 220.2)",
498                "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
499                "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))",
500                "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
501                "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)))",
502                "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))))",
503        };
504
505        for ( i = 0; i < 6; i++ )
506        {
507                geom = lwgeom_over_gserialized(ewkt[i]);
508                CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), LW_TRUE);
509                lwgeom_free(geom);
510        }
511
512        for ( i = 6; i < 12; i++ )
513        {
514                //char *out_ewkt;
515                geom = lwgeom_over_gserialized(ewkt[i]);
516                CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), LW_FALSE);
517                //out_ewkt = lwgeom_to_ewkt(geom, PARSER_CHECK_NONE);
518                //printf("%s\n", out_ewkt);
519                lwgeom_free(geom);
520        }
521
522}
523
524
525static void test_gbox_calculation(void)
526{
527
528        LWGEOM *geom;
529        int i = 0;
530        GBOX *gbox = gbox_new(gflags(0,0,0));
531        BOX3D *box3d;
532
533        char ewkt[][512] =
534        {
535                "POINT(0 0.2)",
536                "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
537                "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))",
538                "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
539                "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)))",
540                "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))))",
541                "POINT(0 220.2)",
542                "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
543                "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))",
544                "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
545                "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)))",
546                "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))))",
547        };
548
549        for ( i = 0; i < 6; i++ )
550        {
551                geom = lwgeom_over_gserialized(ewkt[i]);
552                lwgeom_calculate_gbox(geom, gbox);
553                box3d = lwgeom_compute_box3d(geom);
554                //printf("%g %g\n", gbox->xmin, box3d->xmin);
555                CU_ASSERT_EQUAL(gbox->xmin, box3d->xmin);
556                CU_ASSERT_EQUAL(gbox->xmax, box3d->xmax);
557                CU_ASSERT_EQUAL(gbox->ymin, box3d->ymin);
558                CU_ASSERT_EQUAL(gbox->ymax, box3d->ymax);
559                lwfree(box3d);
560        }
561        lwfree(gbox);
562}
563
564static void test_gserialized_from_lwgeom(void)
565{
566        LWGEOM *geom;
567        GSERIALIZED *g;
568        uint32 type;
569        double *inspect; /* To poke right into the blob. */
570
571        geom = lwgeom_from_ewkt("POINT(0 0.2)", PARSER_CHECK_NONE);
572        g = gserialized_from_lwgeom(geom, 1, 0);
573        type = gserialized_get_type(g);
574        CU_ASSERT_EQUAL( type, POINTTYPE );
575        inspect = (double*)g;
576        CU_ASSERT_EQUAL(inspect[3], 0.2);
577        lwgeom_free(geom);
578        lwfree(g);
579
580        geom = lwgeom_from_ewkt("POLYGON((-1 -1, -1 2.5, 2 2, 2 -1, -1 -1), (0 0, 0 1, 1 1, 1 0, 0 0))", PARSER_CHECK_NONE);
581        g = gserialized_from_lwgeom(geom, 1, 0);
582        type = gserialized_get_type(g);
583        CU_ASSERT_EQUAL( type, POLYGONTYPE );
584        inspect = (double*)g;
585        CU_ASSERT_EQUAL(inspect[9], 2.5);
586        lwgeom_free(geom);
587        lwfree(g);
588
589        geom = lwgeom_from_ewkt("MULTILINESTRING((0 0, 1 1),(0 0.1, 1 1))", PARSER_CHECK_NONE);
590        g = gserialized_from_lwgeom(geom, 1, 0);
591        type = gserialized_get_type(g);
592        CU_ASSERT_EQUAL( type, MULTILINETYPE );
593        inspect = (double*)g;
594        CU_ASSERT_EQUAL(inspect[12], 0.1);
595        lwgeom_free(geom);
596        lwfree(g);
597
598}
599
600static void test_ptarray_point_in_ring(void)
601{
602        LWGEOM *lwg;
603        LWPOLY *poly;
604        POINT2D pt_to_test;
605        POINT2D pt_outside;
606        int result;
607
608        /* Simple containment case */
609        lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE);
610        poly = (LWPOLY*)lwg;
611        pt_to_test.x = 1.05;
612        pt_to_test.y = 1.05;
613        pt_outside.x = 1.2;
614        pt_outside.y = 1.15;
615        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
616        CU_ASSERT_EQUAL(result, LW_TRUE);
617        lwgeom_free(lwg);
618
619        /* Simple noncontainment case */
620        lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE);
621        poly = (LWPOLY*)lwg;
622        pt_to_test.x = 1.05;
623        pt_to_test.y = 1.15;
624        pt_outside.x = 1.2;
625        pt_outside.y = 1.2;
626        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
627        CU_ASSERT_EQUAL(result, LW_FALSE);
628        lwgeom_free(lwg);
629
630        /* Harder noncontainment case */
631        lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE);
632        poly = (LWPOLY*)lwg;
633        pt_to_test.x = 1.05;
634        pt_to_test.y = 0.9;
635        pt_outside.x = 1.2;
636        pt_outside.y = 1.05;
637        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
638        CU_ASSERT_EQUAL(result, LW_FALSE);
639        lwgeom_free(lwg);
640
641        /* Harder containment case */
642        lwg = lwgeom_from_ewkt("POLYGON((0 0, 0 2, 1 2, 0 3, 2 3, 0 4, 3 5, 0 6, 6 10, 6 1, 0 0))", PARSER_CHECK_NONE);
643        poly = (LWPOLY*)lwg;
644        pt_to_test.x = 1.0;
645        pt_to_test.y = 1.0;
646        pt_outside.x = 1.0;
647        pt_outside.y = 10.0;
648        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
649        CU_ASSERT_EQUAL(result, LW_TRUE);
650        lwgeom_free(lwg);
651
652        /* Point on ring at vertex case */
653        lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE);
654        poly = (LWPOLY*)lwg;
655        pt_to_test.x = 1.1;
656        pt_to_test.y = 1.05;
657        pt_outside.x = 1.2;
658        pt_outside.y = 1.05;
659        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
660        CU_ASSERT_EQUAL(result, LW_TRUE);
661        lwgeom_free(lwg);
662
663        /* Point on ring at first vertex case */
664        lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE);
665        poly = (LWPOLY*)lwg;
666        pt_to_test.x = 1.0;
667        pt_to_test.y = 1.0;
668        pt_outside.x = 1.2;
669        pt_outside.y = 1.05;
670        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
671        CU_ASSERT_EQUAL(result, LW_TRUE);
672        lwgeom_free(lwg);
673
674        /* Point on ring between vertexes case */
675        lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE);
676        poly = (LWPOLY*)lwg;
677        pt_to_test.x = 1.0;
678        pt_to_test.y = 1.1;
679        pt_outside.x = 1.2;
680        pt_outside.y = 1.05;
681        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
682        CU_ASSERT_EQUAL(result, LW_TRUE);
683        lwgeom_free(lwg);
684
685        /* Co-linear crossing case for point-in-polygon test, should return LW_TRUE */
686        lwg = lwgeom_from_ewkt("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))", PARSER_CHECK_NONE);
687        poly = (LWPOLY*)lwg;
688        pt_to_test.x = 1.1;
689        pt_to_test.y = 1.05;
690        pt_outside.x = 1.1;
691        pt_outside.y = 1.3;
692        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
693        CU_ASSERT_EQUAL(result, LW_TRUE);
694        lwgeom_free(lwg);
695
696        /* Grazing case for point-in-polygon test, should return LW_FALSE */
697        lwg = lwgeom_from_ewkt("POLYGON((2.0 3.0, 2.0 0.0, 1.0 1.0, 2.0 3.0))", PARSER_CHECK_NONE);
698        lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 2.0, 1.5 1.5, 1.0 1.0))", PARSER_CHECK_NONE);
699        poly = (LWPOLY*)lwg;
700        pt_to_test.x = 1.5;
701        pt_to_test.y = 1.0;
702        pt_outside.x = 1.5;
703        pt_outside.y = 2.0;
704        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
705        CU_ASSERT_EQUAL(result, LW_FALSE);
706        lwgeom_free(lwg);
707
708        /* Grazing case at first point for point-in-polygon test, should return LW_FALSE */
709        lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 2.0 3.0, 2.0 0.0, 1.0 1.0))", PARSER_CHECK_NONE);
710        poly = (LWPOLY*)lwg;
711        pt_to_test.x = 1.0;
712        pt_to_test.y = 0.0;
713        pt_outside.x = 1.0;
714        pt_outside.y = 2.0;
715        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
716        CU_ASSERT_EQUAL(result, LW_FALSE);
717        lwgeom_free(lwg);
718
719        /* Point on vertex of ring */
720        lwg = lwgeom_from_ewkt("POLYGON((-9 50,51 -11,-10 50,-9 50))", PARSER_CHECK_NONE);
721        poly = (LWPOLY*)lwg;
722        pt_to_test.x = -10.0;
723        pt_to_test.y = 50.0;
724        pt_outside.x = -10.2727799838316134;
725        pt_outside.y = -16.9370033133329976;
726        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
727        CU_ASSERT_EQUAL(result, LW_TRUE);
728        lwgeom_free(lwg);
729
730#if 0
731        /* Small polygon and huge distance between outside point and close-but-not-quite-inside point. Should return LW_FALSE. Pretty degenerate case. */
732        lwg = lwgeom_from_ewkt("0103000020E61000000100000025000000ACAD6F91DDB65EC03F84A86D57264540CCABC279DDB65EC0FCE6926B57264540B6DEAA62DDB65EC0A79F6B63572645402E0BE84CDDB65EC065677155572645405D0B1D39DDB65EC0316310425726454082B5DB27DDB65EC060A4E12957264540798BB619DDB65EC0C393A10D57264540D4BC160FDDB65EC0BD0320EE56264540D7AC4E08DDB65EC096C862CC56264540AFD29205DDB65EC02A1F68A956264540363AFA06DDB65EC0722E418656264540B63A780CDDB65EC06E9B0064562645409614E215DDB65EC0E09DA84356264540FF71EF22DDB65EC0B48145265626454036033F33DDB65EC081B8A60C5626454066FB4546DDB65EC08A47A6F7552645409061785BDDB65EC0F05AE0E755264540D4B63772DDB65EC05C86CEDD55264540D2E4C689DDB65EC09B6EBFD95526454082E573A1DDB65EC0C90BD5DB552645401ABE85B8DDB65EC06692FCE35526454039844ECEDDB65EC04D8AF6F155264540928319E2DDB65EC0AD8D570556264540D31055F3DDB65EC02D618F1D56264540343B7A01DEB65EC0EB70CF3956264540920A1A0CDEB65EC03B00515956264540911BE212DEB65EC0E43A0E7B56264540E3F69D15DEB65EC017E4089E562645408D903614DEB65EC0F0D42FC1562645402191B80EDEB65EC0586870E35626454012B84E05DEB65EC09166C80357264540215B41F8DDB65EC08F832B21572645408392F7E7DDB65EC01138C13A57264540F999F0D4DDB65EC0E4A9C14F57264540AC3FB8BFDDB65EC0EED6875F57264540D3DCFEA8DDB65EC04F6C996957264540ACAD6F91DDB65EC03F84A86D57264540", PARSER_CHECK_NONE);
733        poly = (LWPOLY*)lwg;
734        pt_to_test.x = -122.819436560680316;
735        pt_to_test.y = 42.2702301207017328;
736        pt_outside.x = 120.695136159150778;
737        pt_outside.y = 40.6920926049588516;
738        result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
739        CU_ASSERT_EQUAL(result, LW_FALSE);
740        lwgeom_free(lwg);
741#endif
742
743}
744
745static void test_lwpoly_covers_point2d(void)
746{
747        LWPOLY *poly;
748        LWGEOM *lwg;
749        POINT2D pt_to_test;
750        GBOX gbox;
751        int result;
752
753        gbox.flags = gflags(0, 0, 1);
754
755        lwg = lwgeom_from_ewkt("POLYGON((-9 50,51 -11,-10 50,-9 50))", PARSER_CHECK_NONE);
756        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
757        poly = (LWPOLY*)lwg;
758        pt_to_test.x = -10.0;
759        pt_to_test.y = 50.0;
760        result = lwpoly_covers_point2d(poly, &gbox, &pt_to_test);
761        CU_ASSERT_EQUAL(result, LW_TRUE);
762        lwgeom_free(lwg);
763}
764
765
766static void test_lwgeom_distance_sphere(void)
767{
768        LWGEOM *lwg1, *lwg2;
769        GBOX gbox1, gbox2;
770        double d;
771        SPHEROID s;
772
773        /* Init and force spherical */
774        spheroid_init(&s, 6378137.0, 6356752.314245179498);
775        s.a = s.b = s.radius;
776
777        gbox1.flags = gflags(0, 0, 1);
778        gbox2.flags = gflags(0, 0, 1);
779
780        /* Line/line distance, 1 degree apart */
781        lwg1 = lwgeom_from_ewkt("LINESTRING(-30 10, -20 5, -10 3, 0 1)", PARSER_CHECK_NONE);
782        lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", PARSER_CHECK_NONE);
783        lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
784        lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
785        d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
786        CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001);
787        lwgeom_free(lwg1);
788        lwgeom_free(lwg2);
789
790        /* Line/line distance, crossing, 0.0 apart */
791        lwg1 = lwgeom_from_ewkt("LINESTRING(-30 10, -20 5, -10 3, 0 1)", PARSER_CHECK_NONE);
792        lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 20, 5 0, 10 -5)", PARSER_CHECK_NONE);
793        lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
794        lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
795        d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
796        CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
797        lwgeom_free(lwg1);
798        lwgeom_free(lwg2);
799
800        /* Line/point distance, 1 degree apart */
801        lwg1 = lwgeom_from_ewkt("POINT(-4 1)", PARSER_CHECK_NONE);
802        lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", PARSER_CHECK_NONE);
803        lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
804        lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
805        d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
806        CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001);
807        lwgeom_free(lwg1);
808        lwgeom_free(lwg2);
809
810        lwg1 = lwgeom_from_ewkt("POINT(-4 1)", PARSER_CHECK_NONE);
811        lwg2 = lwgeom_from_ewkt("POINT(-4 -1)", PARSER_CHECK_NONE);
812        lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
813        lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
814        d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
815        CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 90.0, 0.00001);
816        lwgeom_free(lwg1);
817        lwgeom_free(lwg2);
818
819        /* Poly/point distance, point inside polygon, 0.0 apart */
820        lwg1 = lwgeom_from_ewkt("POLYGON((-4 1, -3 5, 1 2, 1.5 -5, -4 1))", PARSER_CHECK_NONE);
821        lwg2 = lwgeom_from_ewkt("POINT(-1 -1)", PARSER_CHECK_NONE);
822        lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
823        lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
824        d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
825        CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
826        lwgeom_free(lwg1);
827        lwgeom_free(lwg2);
828
829        /* Poly/point distance, point inside polygon hole, 1 degree apart */
830        lwg1 = lwgeom_from_ewkt("POLYGON((-4 -4, -4 4, 4 4, 4 -4, -4 -4), (-2 -2, -2 2, 2 2, 2 -2, -2 -2))", PARSER_CHECK_NONE);
831        lwg2 = lwgeom_from_ewkt("POINT(-1 -1)", PARSER_CHECK_NONE);
832        lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
833        lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
834        d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
835        CU_ASSERT_DOUBLE_EQUAL(d, 111178.142466, 0.1);
836        lwgeom_free(lwg1);
837        lwgeom_free(lwg2);
838
839        /* Poly/point distance, point on hole boundary, 0.0 apart */
840        lwg1 = lwgeom_from_ewkt("POLYGON((-4 -4, -4 4, 4 4, 4 -4, -4 -4), (-2 -2, -2 2, 2 2, 2 -2, -2 -2))", PARSER_CHECK_NONE);
841        lwg2 = lwgeom_from_ewkt("POINT(2 2)", PARSER_CHECK_NONE);
842        lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
843        lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
844        d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
845        CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
846        lwgeom_free(lwg1);
847        lwgeom_free(lwg2);
848
849        /* Medford test case #1 */
850        lwg1 = lwgeom_from_ewkt("0105000020E610000001000000010200000002000000EF7B8779C7BD5EC0FD20D94B852845400E539C62B9BD5EC0F0A5BE767C284540", PARSER_CHECK_NONE);
851        lwg2 = lwgeom_from_ewkt("0106000020E61000000100000001030000000100000007000000280EC3FB8CCA5EC0A5CDC747233C45402787C8F58CCA5EC0659EA2761E3C45400CED58DF8FCA5EC0C37FAE6E1E3C4540AE97B8E08FCA5EC00346F58B1F3C4540250359FD8ECA5EC05460628E1F3C45403738F4018FCA5EC05DC84042233C4540280EC3FB8CCA5EC0A5CDC747233C4540", PARSER_CHECK_NONE);
852        lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
853        lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
854        d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
855        CU_ASSERT_DOUBLE_EQUAL(d, 23630.8003, 0.1);
856        lwgeom_free(lwg1);
857        lwgeom_free(lwg2);
858
859}
860
861static void test_spheroid_distance(void)
862{
863        GEOGRAPHIC_POINT g1, g2;
864        double d;
865        SPHEROID s;
866
867        /* Init to WGS84 */
868        spheroid_init(&s, 6378137.0, 6356752.314245179498);
869
870        /* One vertical degree */
871        point_set(0.0, 0.0, &g1);
872        point_set(0.0, 1.0, &g2);
873        d = spheroid_distance(&g1, &g2, &s);
874        CU_ASSERT_DOUBLE_EQUAL(d, 110574.388615329, 0.001);
875
876        /* Ten horizontal degrees */
877        point_set(-10.0, 0.0, &g1);
878        point_set(0.0, 0.0, &g2);
879        d = spheroid_distance(&g1, &g2, &s);
880        CU_ASSERT_DOUBLE_EQUAL(d, 1113194.90793274, 0.001);
881
882        /* One horizonal degree */
883        point_set(-1.0, 0.0, &g1);
884        point_set(0.0, 0.0, &g2);
885        d = spheroid_distance(&g1, &g2, &s);
886        CU_ASSERT_DOUBLE_EQUAL(d, 111319.490779, 0.001);
887
888        /* Around world w/ slight bend */
889        point_set(-180.0, 0.0, &g1);
890        point_set(0.0, 1.0, &g2);
891        d = spheroid_distance(&g1, &g2, &s);
892        CU_ASSERT_DOUBLE_EQUAL(d, 19893357.0704483, 0.001);
893
894        /* Up to pole */
895        point_set(-180.0, 0.0, &g1);
896        point_set(0.0, 90.0, &g2);
897        d = spheroid_distance(&g1, &g2, &s);
898        CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7295318, 0.001);
899
900}
901
902static void test_spheroid_area(void)
903{
904        LWGEOM *lwg;
905        GBOX gbox;
906        double a1, a2;
907        SPHEROID s;
908
909        /* Init to WGS84 */
910        spheroid_init(&s, 6378137.0, 6356752.314245179498);
911
912        gbox.flags = gflags(0, 0, 1);
913
914        /* Medford lot test polygon */
915        lwg = lwgeom_from_ewkt("POLYGON((-122.848227067007 42.5007249610493,-122.848309475585 42.5007179884263,-122.848327688675 42.500835880696,-122.848245279942 42.5008428533324,-122.848227067007 42.5007249610493))", PARSER_CHECK_NONE);
916        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
917        a1 = lwgeom_area_sphere(lwg, &gbox, &s);
918        a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
919        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
920        CU_ASSERT_DOUBLE_EQUAL(a1, 89.7211470368, 0.0001); /* sphere */
921        CU_ASSERT_DOUBLE_EQUAL(a2, 89.8684316032, 0.0001); /* spheroid */
922
923        /* Big-ass polygon */
924        lwg = lwgeom_from_ewkt("POLYGON((-2 3, -2 4, -1 4, -1 3, -2 3))", PARSER_CHECK_NONE);
925        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
926        a1 = lwgeom_area_sphere(lwg, &gbox, &s);
927        a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
928        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
929        CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */
930        CU_ASSERT_DOUBLE_EQUAL(a2, 12286574431.9, 10.0); /* spheroid */
931
932        /* One-degree square */
933        lwg = lwgeom_from_ewkt("POLYGON((8.5 2,8.5 1,9.5 1,9.5 2,8.5 2))", PARSER_CHECK_NONE);
934        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
935        a1 = lwgeom_area_sphere(lwg, &gbox, &s);
936        a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
937        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
938        CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */
939        CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
940
941        /* One-degree square *near* dateline */
942        lwg = lwgeom_from_ewkt("POLYGON((179.5 2,179.5 1,178.5 1,178.5 2,179.5 2))", PARSER_CHECK_NONE);
943        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
944        a1 = lwgeom_area_sphere(lwg, &gbox, &s);
945        a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
946        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
947        CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */
948        CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
949
950        /* One-degree square *across* dateline */
951        lwg = lwgeom_from_ewkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", PARSER_CHECK_NONE);
952        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
953        a1 = lwgeom_area_sphere(lwg, &gbox, &s);
954        a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
955        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
956        CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.3679, 10.0); /* sphere */
957        CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
958}
959
960
961/*
962** Used by test harness to register the tests in this file.
963*/
964CU_TestInfo geodetic_tests[] = {
965        PG_TEST(test_signum),
966        PG_TEST(test_gbox_from_spherical_coordinates),
967        PG_TEST(test_gserialized_get_gbox_geocentric),
968        PG_TEST(test_clairaut),
969        PG_TEST(test_edge_intersection),
970        PG_TEST(test_edge_distance_to_point),
971        PG_TEST(test_edge_distance_to_edge),
972        PG_TEST(test_lwgeom_distance_sphere),
973        PG_TEST(test_lwgeom_check_geodetic),
974        PG_TEST(test_gbox_calculation),
975        PG_TEST(test_gserialized_from_lwgeom),
976        PG_TEST(test_spheroid_distance),
977        PG_TEST(test_spheroid_area),
978        PG_TEST(test_lwpoly_covers_point2d),
979        PG_TEST(test_ptarray_point_in_ring),
980        CU_TEST_INFO_NULL
981};
982CU_SuiteInfo geodetic_suite = {"Geodetic Suite",  NULL,  NULL, geodetic_tests};
Note: See TracBrowser for help on using the repository browser.