capi/geostest.c

   1 /************************************************************************
   2  *
   3  * $Id$
   4  *
   5  * Test for C-Wrapper of GEOS library
   6  *
   7  * Copyright (C) 2005 Refractions Research Inc.
   8  *
   9  * This is free software; you can redistribute and/or modify it under
  10  * the terms of the GNU Lesser General Public Licence as published
  11  * by the Free Software Foundation.
  12  * See the COPYING file for more information.
  13  *
  14  * Author: Sandro Santilli <strk@refractions.net>
  15  *
  16  ***********************************************************************/
  17
  18 #define _GNU_SOURCE
  19
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22 #include <stdarg.h>
  23
  24 #include "geos_c.h"
  25
  26 #define MAXWKTLEN 1047551
  27
  28 void
  29 usage(char *me)
  30 {
  31         fprintf(stderr, "Usage: %s <wktfile>\n", me);
  32         exit(1);
  33 }
  34
  35 void
  36 notice(const char *fmt, ...) {
  37         va_list ap;
  38
  39         fprintf( stdout, "NOTICE: ");
  40
  41         va_start (ap, fmt);
  42         vfprintf( stdout, fmt, ap);
  43         va_end(ap);
  44         fprintf( stdout, "\n" );
  45 }
  46
  47 void
  48 log_and_exit(const char *fmt, ...) {
  49         va_list ap;
  50
  51         fprintf( stdout, "ERROR: ");
  52
  53         va_start (ap, fmt);
  54         vfprintf( stdout, fmt, ap);
  55         va_end(ap);
  56         fprintf( stdout, "\n" );
  57         exit(1);
  58 }
  59
  60 GEOSGeometry*
  61 fineGrainedReconstructionTest(const GEOSGeometry* g1)
  62 {
  63         GEOSCoordSequence* cs;
  64         GEOSGeometry* g2;
  65         GEOSGeometry* shell;
  66         const GEOSGeometry* gtmp;
  67         GEOSGeometry**geoms;
  68         unsigned int ngeoms, i;
  69         int type;
  70
  71         /* Geometry reconstruction from CoordSeq */
  72         type = GEOSGeomTypeId(g1);
  73         switch ( type )
  74         {
  75                 case GEOS_POINT:
  76                         cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(g1));
  77                         g2 = GEOSGeom_createPoint(cs);
  78                         return g2;
  79                         break;
  80                 case GEOS_LINESTRING:
  81                         cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(g1));
  82                         g2 = GEOSGeom_createLineString(cs);
  83                         return g2;
  84                         break;
  85                 case GEOS_LINEARRING:
  86                         cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(g1));
  87                         g2 = GEOSGeom_createLinearRing(cs);
  88                         return g2;
  89                         break;
  90                 case GEOS_POLYGON:
  91                         gtmp = GEOSGetExteriorRing(g1);
  92                         cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(gtmp));
  93                         shell = GEOSGeom_createLinearRing(cs);
  94                         ngeoms = GEOSGetNumInteriorRings(g1);
  95                         geoms = malloc(ngeoms*sizeof(GEOSGeometry*));
  96                         for (i=0; i<ngeoms; i++)
  97                         {
  98                                 gtmp = GEOSGetInteriorRingN(g1, i);
  99                                 cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(gtmp));
 100                                 geoms[i] = GEOSGeom_createLinearRing(cs);
 101                         }
 102                         g2 = GEOSGeom_createPolygon(shell, geoms, ngeoms);
 103                         free(geoms);
 104                         return g2;
 105                         break;
 106                 case GEOS_MULTIPOINT:
 107                 case GEOS_MULTILINESTRING:
 108                 case GEOS_MULTIPOLYGON:
 109                 case GEOS_GEOMETRYCOLLECTION:
 110                         ngeoms = GEOSGetNumGeometries(g1);
 111                         geoms = malloc(ngeoms*sizeof(GEOSGeometry*));
 112                         for (i=0; i<ngeoms; i++)
 113                         {
 114                                 gtmp = GEOSGetGeometryN(g1, i);
 115                                 geoms[i] = fineGrainedReconstructionTest(gtmp);
 116                         }
 117                         g2 = GEOSGeom_createCollection(type, geoms, ngeoms);
 118                         free(geoms);
 119                         return g2;
 120                         break;
 121                 default:
 122                         log_and_exit("Unknown geometry type %d\n", type);
 123                         return NULL;
 124         }
 125 }
 126
 127 void
 128 printHEX(FILE *where, const unsigned char *bytes, size_t n)
 129 {
 130         static char hex[] = "0123456789ABCDEF";
 131         int i;
 132
 133         for (i=0; i<n; i++)
 134         {
 135                 fprintf(where, "%c%c", hex[bytes[i]>>4], hex[bytes[i]&0x0F]);
 136         }
 137 }
 138
 139 int
 140 do_all(char *inputfile)
 141 {
 142         GEOSGeometry* g1;
 143         GEOSGeometry* g2;
 144         GEOSGeometry* g3;
 145         GEOSGeometry* g4;
 146         const GEOSGeometry **gg;
 147         unsigned int npoints, ndims;
 148         static char wkt[MAXWKTLEN];
 149         FILE *input;
 150         char *ptr;
 151         unsigned char* uptr;
 152         size_t size;
 153         double dist, area;
 154
 155         input = fopen(inputfile, "r");
 156         if ( ! input ) { perror("fopen"); exit(1); }
 157
 158         size = fread(wkt, 1, MAXWKTLEN-1, input);
 159         fclose(input);
 160         if ( ! size ) { perror("fread"); exit(1); }
 161         if ( size == MAXWKTLEN-1 ) { perror("WKT input too big!"); exit(1); }
 162         wkt[size] = '\0'; /* ensure it is null terminated */
 163
 164         /* WKT input */
 165         g1 = GEOSGeomFromWKT(wkt);
 166
 167         /* WKT output */
 168         ptr = GEOSGeomToWKT(g1);
 169         printf("Input (WKT): %s\n", ptr);
 170         free(ptr);
 171
 172         /* WKB output */
 173         uptr = GEOSGeomToWKB_buf(g1, &size);
 174         printf("Input (WKB): "); printHEX(stdout, uptr, size); putchar('\n');
 175
 176         /* WKB input */
 177         g2 = GEOSGeomFromWKB_buf(uptr, size); free(uptr);
 178         if ( ! GEOSEquals(g1, g2) ) log_and_exit("Round WKB conversion failed");
 179         GEOSGeom_destroy(g2);
 180
 181         /* Size and dimension */
 182         npoints = GEOSGetNumCoordinates(g1);
 183         ndims = GEOSGeom_getDimensions(g1);
 184         printf("Geometry coordinates: %dx%d\n", npoints, ndims);
 185
 186         /* Geometry fine-grained deconstruction/reconstruction  test */
 187         g2 = fineGrainedReconstructionTest(g1);
 188         if ( ! GEOSEquals(g1, g2) )
 189         {
 190                 log_and_exit("Reconstruction test failed\n");
 191         }
 192         GEOSGeom_destroy(g2);
 193
 194         /* Unary predicates */
 195         if ( GEOSisEmpty(g1) ) printf("isEmpty\n");
 196         if ( GEOSisValid(g1) ) printf("isValid\n");
 197         if ( GEOSisSimple(g1) ) printf("isSimple\n");
 198         if ( GEOSisRing(g1) ) printf("isRing\n");
 199
 200         /* Convex Hull */
 201         g2 = GEOSConvexHull(g1);
 202         if ( ! g2 )
 203         {
 204                 log_and_exit("GEOSConvexHull() raised an exception");
 205         }
 206         ptr = GEOSGeomToWKT(g2);
 207         printf("ConvexHull: %s\n", ptr);
 208         free(ptr);
 209
 210         /* Buffer */
 211         GEOSGeom_destroy(g1);
 212         g1 = GEOSBuffer(g2, 100, 30);
 213         if ( ! g1 )
 214         {
 215                 log_and_exit("GEOSBuffer() raised an exception");
 216         }
 217         ptr = GEOSGeomToWKT(g1);
 218         printf("Buffer: %s\n", ptr);
 219         free(ptr);
 220
 221
 222         /* Intersection */
 223         g3 = GEOSIntersection(g1, g2);
 224         if ( ! GEOSEquals(g3, g2) )
 225         {
 226                 GEOSGeom_destroy(g1);
 227                 GEOSGeom_destroy(g2);
 228                 GEOSGeom_destroy(g3);
 229                 log_and_exit("Intersection(g, Buffer(g)) didn't return g");
 230         }
 231         ptr = GEOSGeomToWKT(g3);
 232         printf("Intersection: %s\n", ptr);
 233         GEOSGeom_destroy(g3);
 234         free(ptr);
 235
 236         /* Difference */
 237         g3 = GEOSDifference(g1, g2);
 238         ptr = GEOSGeomToWKT(g3);
 239         printf("Difference: %s\n", ptr);
 240         GEOSGeom_destroy(g3);
 241         free(ptr);
 242
 243         /* SymDifference */
 244         g3 = GEOSSymDifference(g1, g2);
 245         ptr = GEOSGeomToWKT(g3);
 246         printf("SymDifference: %s\n", ptr);
 247         free(ptr);
 248
 249         /* Boundary */
 250         g4 = GEOSBoundary(g3);
 251         ptr = GEOSGeomToWKT(g4);
 252         printf("Boundary: %s\n", ptr);
 253         GEOSGeom_destroy(g3);
 254         GEOSGeom_destroy(g4);
 255         free(ptr);
 256
 257         /* Union */
 258         g3 = GEOSUnion(g1, g2);
 259         if ( ! GEOSEquals(g3, g1) )
 260         {
 261                 GEOSGeom_destroy(g1);
 262                 GEOSGeom_destroy(g2);
 263                 GEOSGeom_destroy(g3);
 264                 log_and_exit("Union(g, Buffer(g)) didn't return Buffer(g)");
 265         }
 266         ptr = GEOSGeomToWKT(g3);
 267         printf("Union: %s\n", ptr);
 268         free(ptr);
 269
 270         /* PointOnSurcace */
 271         g4 = GEOSPointOnSurface(g3);
 272         ptr = GEOSGeomToWKT(g4);
 273         printf("PointOnSurface: %s\n", ptr);
 274         GEOSGeom_destroy(g3);
 275         GEOSGeom_destroy(g4);
 276         free(ptr);
 277
 278         /* Centroid */
 279         g3 = GEOSGetCentroid(g2);
 280         ptr = GEOSGeomToWKT(g3);
 281         printf("Centroid: %s\n", ptr);
 282         GEOSGeom_destroy(g3);
 283         free(ptr);
 284
 285         /* Relate (and RelatePattern )*/
 286         ptr = GEOSRelate(g1, g2);
 287         if ( ! GEOSRelatePattern(g1, g2, ptr) )
 288         {
 289                 GEOSGeom_destroy(g1);
 290                 GEOSGeom_destroy(g2);
 291                 free(ptr);
 292                 log_and_exit("! RelatePattern(g1, g2, Relate(g1, g2))");
 293         }
 294         printf("Relate: %s\n", ptr);
 295         free(ptr);
 296
 297         /* Polygonize */
 298         gg = (const GEOSGeometry**)malloc(2*sizeof(GEOSGeometry*));
 299         gg[0] = g1;
 300         gg[1] = g2;
 301         g3 = GEOSPolygonize(gg, 2);
 302         free(gg);
 303         if ( ! g3 )
 304         {
 305                 log_and_exit("Exception running GEOSPolygonize");
 306         }
 307         ptr = GEOSGeomToWKT(g3);
 308         GEOSGeom_destroy(g3);
 309         printf("Polygonize: %s\n", ptr);
 310         free(ptr);
 311
 312         /* LineMerge */
 313         g3 = GEOSLineMerge(g1);
 314         if ( ! g3 )
 315         {
 316                 log_and_exit("Exception running GEOSLineMerge");
 317         }
 318         ptr = GEOSGeomToWKT(g3);
 319         printf("LineMerge: %s\n", ptr);
 320         free(ptr);
 321         GEOSGeom_destroy(g3);
 322
 323         /* Binary predicates */
 324         if ( GEOSIntersects(g1, g2) ) printf("Intersect\n");
 325         if ( GEOSDisjoint(g1, g2) ) printf("Disjoint\n");
 326         if ( GEOSTouches(g1, g2) ) printf("Touches\n");
 327         if ( GEOSCrosses(g1, g2) ) printf("Crosses\n");
 328         if ( GEOSWithin(g1, g2) ) printf("Within\n");
 329         if ( GEOSContains(g1, g2) ) printf("Contains\n");
 330         if ( GEOSOverlaps(g1, g2) ) printf("Overlaps\n");
 331
 332         /* Distance */
 333         if ( GEOSDistance(g1, g2, &dist) ) printf("Distance: %g\n", dist);
 334
 335         /* Area */
 336         if ( GEOSArea(g1, &area) ) printf("Area 1: %g\n", area);
 337         if ( GEOSArea(g2, &area) ) printf("Area 2: %g\n", area);
 338
 339         GEOSGeom_destroy(g2);
 340
 341         /* Simplify */
 342         g3 = GEOSSimplify(g1, 0.5);
 343         ptr = GEOSGeomToWKT(g3);
 344         printf("Simplify: %s\n", ptr);
 345         free(ptr);
 346         GEOSGeom_destroy(g3);
 347
 348         /* Topology Preserve Simplify */
 349         g3 = GEOSTopologyPreserveSimplify(g1, 0.5);
 350         ptr = GEOSGeomToWKT(g3);
 351         printf("Simplify: %s\n", ptr);
 352         free(ptr);
 353         GEOSGeom_destroy(g3);
 354
 355         GEOSGeom_destroy(g1);
 356
 357         return 0;
 358 }
 359
 360 int
 361 main(int argc, char **argv)
 362 {
 363         int i, n=1;
 364
 365         initGEOS(notice, log_and_exit);
 366         printf("GEOS version %s\n", GEOSversion());
 367
 368         if ( argc < 2 ) usage(argv[0]);
 369         if ( argc > 2 ) n=atoi(argv[2]);
 370         if ( ! n ) n=1;
 371
 372         for (i=0; i<n; i++) {
 373                 putc('.', stderr); fflush(stderr);
 374                 do_all(argv[1]);
 375                 putc('+', stderr); fflush(stderr);
 376         }
 377         putc('\n', stderr);
 378
 379         finishGEOS();
 380
 381         return EXIT_SUCCESS;
 382 }
 383