1 /**********************************************************************
2 *
3 * GEOS - Geometry Engine Open Source
4 * http://geos.osgeo.org
5 *
6 * Copyright (C) 2001-2002 Vivid Solutions Inc.
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 *********************************************************************
15 *
16 * This file should document by example usage of the GEOS library.
17 * It could actually be a live discuss-by-example board for
18 * architectural design choices.
19 *
20 * --strk;
21 *
22 * DEBUGGING TIPS:
23 * use -D__USE_MALLOC at compile time for gcc 2.91, 2.95, 3.0 and 3.1
24 * and GLIBCXX_FORCE_NEW or GLIBCPP_FORCE_NEW at run time with gcc 3.2.2+
25 * to force libstdc++ avoid caching memory. This should remove some
26 * obscure reports from memory checkers like valgrind.
27 *
28 **********************************************************************/
29
30 #include <geos/geom/PrecisionModel.h>
31 #include <geos/geom/GeometryFactory.h>
32 #include <geos/geom/Geometry.h>
33 #include <geos/geom/Point.h>
34 #include <geos/geom/LinearRing.h>
35 #include <geos/geom/LineString.h>
36 #include <geos/geom/Polygon.h>
37 #include <geos/geom/GeometryCollection.h>
38 #include <geos/geom/Coordinate.h>
39 #include <geos/geom/CoordinateSequence.h>
40 #include <geos/geom/CoordinateArraySequence.h>
41 #include <geos/geom/IntersectionMatrix.h>
42 #include <geos/io/WKBReader.h>
43 #include <geos/io/WKBWriter.h>
44 #include <geos/io/WKTWriter.h>
45 #include <geos/util/GeometricShapeFactory.h>
46 #include <geos/geom/util/SineStarFactory.h>
47 #include <geos/util/GEOSException.h>
48 #include <geos/util/IllegalArgumentException.h>
49 #include <geos/operation/linemerge/LineMerger.h>
50 #include <geos/operation/polygonize/Polygonizer.h>
51 #include <geos/constants.h>
52 #include <vector>
53 #include <sstream>
54 #include <iomanip>
55 #include <cstdlib> // exit()
56
57 // Set to 0 to skip section
58 #define GEOMETRIC_SHAPES 1
59 #define RELATIONAL_OPERATORS 1
60 #define COMBINATIONS 1
61 #define UNARY_OPERATIONS 1
62 #define LINEMERGE 1
63 #define POLYGONIZE 1
64
65
66 using namespace std;
67 using namespace geos;
68 using namespace geos::geom;
69 using namespace geos::operation::polygonize;
70 using namespace geos::operation::linemerge;
71 using geos::util::GEOSException;
72 using geos::util::IllegalArgumentException;
73
74
75 // Prototypes
76 void wkt_print_geoms(vector<const Geometry*>* geoms);
77
78
79 // This object will be used to construct our geometries.
80 // It might be bypassed by directly call geometry constructors,
81 // but that would be boring because you'd need to specify
82 // a PrecisionModel and a SRID everytime: those infos are
83 // cached inside a GeometryFactory object.
84 GeometryFactory::Ptr global_factory;
85
86 //#define DEBUG_STREAM_STATE 1
87
88
89 //
90 // This function tests writing and reading WKB
91 // TODO:
92 // - compare input and output geometries for equality
93 // - remove debugging lines (on stream state)
94 //
95 void
WKBtest(vector<const Geometry * > * geoms)96 WKBtest(vector<const Geometry*>* geoms)
97 {
98 stringstream s(ios_base::binary | ios_base::in | ios_base::out);
99 io::WKBReader wkbReader(*global_factory);
100 io::WKBWriter wkbWriter;
101 Geometry* gout;
102
103 #if DEBUG_STREAM_STATE
104 cout << "WKBtest: machine byte order: " << BYTE_ORDER << endl;
105 #endif
106
107
108 size_t ngeoms = geoms->size();
109 for(unsigned int i = 0; i < ngeoms; ++i) {
110 const Geometry* gin = (*geoms)[i];
111
112 #if DEBUG_STREAM_STATE
113 cout << "State of stream before WRITE: ";
114 cout << "p:" << s.tellp() << " g:" << s.tellg() <<
115 " good:" << s.good() <<
116 " eof:" << s.eof() <<
117 " bad:" << s.bad() <<
118 " fail:" << s.fail() << endl;
119 #endif
120
121 #if DEBUG_STREAM_STATE
122 cout << "State of stream after SEEKP: ";
123 cout << "p:" << s.tellp() << " g:" << s.tellg() <<
124 " good:" << s.good() <<
125 " eof:" << s.eof() <<
126 " bad:" << s.bad() <<
127 " fail:" << s.fail() << endl;
128 #endif
129
130 wkbWriter.write(*gin, s);
131 #if DEBUG_STREAM_STATE
132 cout << "wkbWriter wrote and reached ";
133 cout << "p:" << s.tellp() << " g:" << s.tellg() << endl;
134
135 cout << "State of stream before DUMP: ";
136 cout << "p:" << s.tellp() << " g:" << s.tellg() <<
137 " good:" << s.good() <<
138 " eof:" << s.eof() <<
139 " bad:" << s.bad() <<
140 " fail:" << s.fail() << endl;
141 #endif
142
143 #if DEBUG_STREAM_STATE
144 cout << "State of stream after DUMP: ";
145 cout << "p:" << s.tellp() << " g:" << s.tellg() <<
146 " good:" << s.good() <<
147 " eof:" << s.eof() <<
148 " bad:" << s.bad() <<
149 " fail:" << s.fail() << endl;
150 #endif
151
152 s.seekg(0, ios::beg); // rewind reader pointer
153
154 #if DEBUG_STREAM_STATE
155 cout << "State of stream before READ: ";
156 cout << "p:" << s.tellp() << " g:" << s.tellg() <<
157 " good:" << s.good() <<
158 " eof:" << s.eof() <<
159 " bad:" << s.bad() <<
160 " fail:" << s.fail() << endl;
161 #endif
162
163 gout = wkbReader.read(s).release();
164
165 #if DEBUG_STREAM_STATE
166 cout << "State of stream after READ: ";
167 cout << "p:" << s.tellp() << " g:" << s.tellg() <<
168 " good:" << s.good() <<
169 " eof:" << s.eof() <<
170 " bad:" << s.bad() <<
171 " fail:" << s.fail() << endl;
172 #endif
173
174 const_cast<Geometry*>(gin)->normalize();
175 gout->normalize();
176 int failed = gin->compareTo(gout);
177 if(failed) {
178 cout << "{" << i << "} (WKB) ";
179 }
180 else {
181 cout << "[" << i << "] (WKB) ";
182 }
183
184 io::WKBReader::printHEX(s, cout);
185 cout << endl;
186
187 if(failed) {
188 io::WKTWriter wkt;
189 cout << " IN: " << wkt.write(gin) << endl;
190 cout << " OUT: " << wkt.write(gout) << endl;
191 }
192
193 s.seekp(0, ios::beg); // rewind writer pointer
194
195 delete gout;
196 }
197
198 }
199
200
201 // This function will print given geometries in WKT
202 // format to stdout. As a side-effect, will test WKB
203 // output and input, using the WKBtest function.
204 void
wkt_print_geoms(vector<const Geometry * > * geoms)205 wkt_print_geoms(vector<const Geometry*>* geoms)
206 {
207 WKBtest(geoms); // test WKB parser
208
209 // WKT-print given geometries
210 io::WKTWriter* wkt = new io::WKTWriter();
211 for(unsigned int i = 0; i < geoms->size(); i++) {
212 const Geometry* g = (*geoms)[i];
213 string tmp = wkt->write(g);
214 cout << "[" << i << "] (WKT) " << tmp << endl;
215 }
216 delete wkt;
217 }
218
219 // This is the simpler geometry you can get: a point.
220 Point*
create_point(double x,double y)221 create_point(double x, double y)
222 {
223 Coordinate c(x, y);
224 Point* p = global_factory->createPoint(c);
225 return p;
226 }
227
228 // This function will create a LinearString
229 // geometry with the shape of the letter U
230 // having top-left corner at given coordinates
231 // and 'side' height and width
232 LineString*
create_ushaped_linestring(double xoffset,double yoffset,double side)233 create_ushaped_linestring(double xoffset, double yoffset, double side)
234 {
235 // We will use a coordinate list to build the linestring
236 CoordinateArraySequence* cl = new CoordinateArraySequence();
237
238 cl->add(Coordinate(xoffset, yoffset));
239 cl->add(Coordinate(xoffset, yoffset + side));
240 cl->add(Coordinate(xoffset + side, yoffset + side));
241 cl->add(Coordinate(xoffset + side, yoffset));
242
243 // Now that we have a CoordinateSequence we can create
244 // the linestring.
245 // The newly created LineString will take ownership
246 // of the CoordinateSequence.
247 LineString* ls = global_factory->createLineString(cl);
248
249 // This is what you do if you want the new LineString
250 // to make a copy of your CoordinateSequence:
251 // LineString *ls = global_factory->createLineString(*cl);
252
253 return ls; // our LineString
254 }
255
256 // This function will create a LinearRing
257 // geometry rapresenting a square with the given origin
258 // and side
259 LinearRing*
create_square_linearring(double xoffset,double yoffset,double side)260 create_square_linearring(double xoffset, double yoffset, double side)
261 {
262 // We will use a coordinate list to build the linearring
263 CoordinateArraySequence* cl = new CoordinateArraySequence();
264
265 cl->add(Coordinate(xoffset, yoffset));
266 cl->add(Coordinate(xoffset, yoffset + side));
267 cl->add(Coordinate(xoffset + side, yoffset + side));
268 cl->add(Coordinate(xoffset + side, yoffset));
269 cl->add(Coordinate(xoffset, yoffset));
270
271 // Now that we have a CoordinateSequence we can create
272 // the linearring.
273 // The newly created LinearRing will take ownership
274 // of the CoordinateSequence.
275 LinearRing* lr = global_factory->createLinearRing(cl);
276
277 // This is what you do if you want the new LinearRing
278 // to make a copy of your CoordinateSequence:
279 // LinearRing *lr = global_factory->createLinearRing(*cl);
280
281 return lr; // our LinearRing
282 }
283
284 // This function will create a Polygon
285 // geometry rapresenting a square with the given origin
286 // and side and with a central hole 1/3 sided.
287 Polygon*
create_square_polygon(double xoffset,double yoffset,double side)288 create_square_polygon(double xoffset, double yoffset, double side)
289 {
290 // We need a LinearRing for the polygon shell
291 LinearRing* outer = create_square_linearring(xoffset, yoffset, side);
292
293 // And another for the hole
294 LinearRing* inner = create_square_linearring(xoffset + (side / 3),
295 yoffset + (side / 3), (side / 3));
296
297 // If we need to specify any hole, we do it using
298 // a vector of LinearRing pointers
299 vector<LinearRing*>* holes = new vector<LinearRing*>;
300
301 // We add the newly created geometry to the vector
302 // of holes.
303 holes->push_back(inner);
304
305 // And finally we call the polygon constructor.
306 // Both the outer LinearRing and the vector of holes
307 // will be referenced by the resulting Polygon object,
308 // thus we CANNOT delete them, neither the holes, nor
309 // the vector containing their pointers, nor the outer
310 // LinearRing. Everything will be deleted at Polygon
311 // deletion time (this is inconsistent with LinearRing
312 // behaviour... what should we do?).
313 Polygon* poly = global_factory->createPolygon(outer, holes);
314
315 return poly;
316 }
317
318 //
319 // This function will create a GeometryCollection
320 // containing copies of all Geometries in given vector.
321 //
322 GeometryCollection*
create_simple_collection(vector<const Geometry * > * geoms)323 create_simple_collection(vector<const Geometry*>* geoms)
324 {
325 return global_factory->createGeometryCollection(*geoms);
326 // if you wanted to transfer ownership of vector end
327 // its elements you should have call:
328 // return global_factory->createGeometryCollection(geoms);
329 }
330
331 //
332 // This function uses GeometricShapeFactory to render
333 // a circle having given center and radius
334 //
335 Polygon*
create_circle(double centerX,double centerY,double radius)336 create_circle(double centerX, double centerY, double radius)
337 {
338 geos::util::GeometricShapeFactory shapefactory(global_factory.get());
339 shapefactory.setCentre(Coordinate(centerX, centerY));
340 shapefactory.setSize(radius);
341 // same as:
342 // shapefactory.setHeight(radius);
343 // shapefactory.setWidth(radius);
344 return shapefactory.createCircle().release();
345 }
346
347 //
348 // This function uses GeometricShapeFactory to render
349 // an ellipse having given center and axis size
350 //
351 Polygon*
create_ellipse(double centerX,double centerY,double width,double height)352 create_ellipse(double centerX, double centerY, double width, double height)
353 {
354 geos::util::GeometricShapeFactory shapefactory(global_factory.get());
355 shapefactory.setCentre(Coordinate(centerX, centerY));
356 shapefactory.setHeight(height);
357 shapefactory.setWidth(width);
358 return shapefactory.createCircle().release();
359 }
360
361 //
362 // This function uses GeometricShapeFactory to render
363 // a rectangle having lower-left corner at given coordinates
364 // and given sizes.
365 //
366 Polygon*
create_rectangle(double llX,double llY,double width,double height)367 create_rectangle(double llX, double llY, double width, double height)
368 {
369 geos::util::GeometricShapeFactory shapefactory(global_factory.get());
370 shapefactory.setBase(Coordinate(llX, llY));
371 shapefactory.setHeight(height);
372 shapefactory.setWidth(width);
373 shapefactory.setNumPoints(4); // we don't need more then 4 points for a rectangle...
374 // can use setSize for a square
375 return shapefactory.createRectangle().release();
376 }
377
378 //
379 // This function uses GeometricShapeFactory to render
380 // an arc having lower-left corner at given coordinates,
381 // given sizes and given angles.
382 //
383 LineString*
create_arc(double llX,double llY,double width,double height,double startang,double endang)384 create_arc(double llX, double llY, double width, double height, double startang, double endang)
385 {
386 geos::util::GeometricShapeFactory shapefactory(global_factory.get());
387 shapefactory.setBase(Coordinate(llX, llY));
388 shapefactory.setHeight(height);
389 shapefactory.setWidth(width);
390 // shapefactory.setNumPoints(100); // the default (100 pts)
391 // can use setSize for a square
392 return shapefactory.createArc(startang, endang).release();
393 }
394
395 unique_ptr<Polygon>
create_sinestar(double cx,double cy,double size,int nArms,double armLenRat)396 create_sinestar(double cx, double cy, double size, int nArms, double armLenRat)
397 {
398 geos::geom::util::SineStarFactory fact(global_factory.get());
399 fact.setCentre(Coordinate(cx, cy));
400 fact.setSize(size);
401 fact.setNumPoints(nArms * 5);
402 fact.setArmLengthRatio(armLenRat);
403 fact.setNumArms(nArms);
404 return fact.createSineStar();
405 }
406
407 // Start reading here
408 void
do_all()409 do_all()
410 {
411 vector<const Geometry*>* geoms = new vector<const Geometry*>;
412 vector<const Geometry*>* newgeoms;
413
414 // Define a precision model using 0,0 as the reference origin
415 // and 2.0 as coordinates scale.
416 PrecisionModel* pm = new PrecisionModel(2.0, 0, 0);
417
418 // Initialize global factory with defined PrecisionModel
419 // and a SRID of -1 (undefined).
420 global_factory = GeometryFactory::create(pm, -1);
421
422 // We do not need PrecisionMode object anymore, it has
423 // been copied to global_factory private storage
424 delete pm;
425
426 ////////////////////////////////////////////////////////////////////////
427 // GEOMETRY CREATION
428 ////////////////////////////////////////////////////////////////////////
429
430 // Read function bodies to see the magic behind them
431 geoms->push_back(create_point(150, 350));
432 geoms->push_back(create_square_linearring(0, 0, 100));
433 geoms->push_back(create_ushaped_linestring(60, 60, 100));
434 geoms->push_back(create_square_linearring(0, 0, 100));
435 geoms->push_back(create_square_polygon(0, 200, 300));
436 geoms->push_back(create_square_polygon(0, 250, 300));
437 geoms->push_back(create_simple_collection(geoms));
438
439 #if GEOMETRIC_SHAPES
440 // These ones use a GeometricShapeFactory
441 geoms->push_back(create_circle(0, 0, 10));
442 geoms->push_back(create_ellipse(0, 0, 8, 12));
443 geoms->push_back(create_rectangle(-5, -5, 10, 10)); // a square
444 geoms->push_back(create_rectangle(-5, -5, 10, 20)); // a rectangle
445 // The upper-right quarter of a vertical ellipse
446 geoms->push_back(create_arc(0, 0, 10, 20, 0, MATH_PI / 2));
447 geoms->push_back(create_sinestar(10, 10, 100, 5, 2).release()); // a sine star
448 #endif
449
450 // Print all geoms.
451 cout << "--------HERE ARE THE BASE GEOMS ----------" << endl;
452 wkt_print_geoms(geoms);
453
454
455 #if UNARY_OPERATIONS
456
457 ////////////////////////////////////////////////////////////////////////
458 // UNARY OPERATIONS
459 ////////////////////////////////////////////////////////////////////////
460
461 /////////////////////////////////////////////
462 // CENTROID
463 /////////////////////////////////////////////
464
465 // Find centroid of each base geometry
466 newgeoms = new vector<const Geometry*>;
467 for(unsigned int i = 0; i < geoms->size(); i++) {
468 const Geometry* g = (*geoms)[i];
469 newgeoms->push_back(g->getCentroid().release());
470 }
471
472 // Print all convex hulls
473 cout << endl << "------- AND HERE ARE THEIR CENTROIDS -----" << endl;
474 wkt_print_geoms(newgeoms);
475
476 // Delete the centroids
477 for(unsigned int i = 0; i < newgeoms->size(); i++) {
478 delete(*newgeoms)[i];
479 }
480 delete newgeoms;
481
482 /////////////////////////////////////////////
483 // BUFFER
484 /////////////////////////////////////////////
485
486 newgeoms = new vector<const Geometry*>;
487 for(unsigned int i = 0; i < geoms->size(); i++) {
488 const Geometry* g = (*geoms)[i];
489 try {
490 Geometry* g2 = g->buffer(10).release();
491 newgeoms->push_back(g2);
492 }
493 catch(const GEOSException& exc) {
494 cerr << "GEOS Exception: geometry " << i << "->buffer(10): " << exc.what() << "\n";
495 }
496 }
497
498 cout << endl << "--------HERE COMES THE BUFFERED GEOMS ----------" << endl;
499 wkt_print_geoms(newgeoms);
500
501 for(unsigned int i = 0; i < newgeoms->size(); i++) {
502 delete(*newgeoms)[i];
503 }
504 delete newgeoms;
505
506 /////////////////////////////////////////////
507 // CONVEX HULL
508 /////////////////////////////////////////////
509
510 // Make convex hulls of geometries
511 newgeoms = new vector<const Geometry*>;
512 for(unsigned int i = 0; i < geoms->size(); i++) {
513 const Geometry* g = (*geoms)[i];
514 newgeoms->push_back(g->convexHull().release());
515 }
516
517 // Print all convex hulls
518 cout << endl << "--------HERE COMES THE HULLS----------" << endl;
519 wkt_print_geoms(newgeoms);
520
521 // Delete the hulls
522 for(unsigned int i = 0; i < newgeoms->size(); i++) {
523 delete(*newgeoms)[i];
524 }
525 delete newgeoms;
526
527 #endif // UNARY_OPERATIONS
528
529 #if RELATIONAL_OPERATORS
530
531 ////////////////////////////////////////////////////////////////////////
532 // RELATIONAL OPERATORS
533 ////////////////////////////////////////////////////////////////////////
534
535 cout << "-------------------------------------------------------------------------------" << endl;
536 cout << "RELATIONAL OPERATORS" << endl;
537 cout << "-------------------------------------------------------------------------------" << endl;
538
539 /////////////////////////////////////////////
540 // DISJOINT
541 /////////////////////////////////////////////
542
543 cout << endl;
544 cout << " DISJOINT ";
545 for(unsigned int i = 0; i < geoms->size(); i++) {
546 cout << "\t[" << i << "]";
547 }
548 cout << endl;
549 for(unsigned int i = 0; i < geoms->size(); i++) {
550 const Geometry* g1 = (*geoms)[i];
551 cout << " [" << i << "]\t";
552 for(unsigned int j = 0; j < geoms->size(); j++) {
553 const Geometry* g2 = (*geoms)[j];
554 try {
555 if(g1->disjoint(g2)) {
556 cout << " 1\t";
557 }
558 else {
559 cout << " 0\t";
560 }
561 }
562 // Geometry Collection is not a valid argument
563 catch(const IllegalArgumentException&) {
564 cout << " X\t";
565 }
566 catch(const std::exception& exc) {
567 cerr << exc.what() << endl;
568 }
569 }
570 cout << endl;
571 }
572
573 /////////////////////////////////////////////
574 // TOUCHES
575 /////////////////////////////////////////////
576
577 cout << endl;
578 cout << " TOUCHES ";
579 for(unsigned int i = 0; i < geoms->size(); i++) {
580 cout << "\t[" << i << "]";
581 }
582 cout << endl;
583 for(unsigned int i = 0; i < geoms->size(); i++) {
584 const Geometry* g1 = (*geoms)[i];
585 cout << " [" << i << "]\t";
586 for(unsigned int j = 0; j < geoms->size(); j++) {
587 const Geometry* g2 = (*geoms)[j];
588 try {
589 if(g1->touches(g2)) {
590 cout << " 1\t";
591 }
592 else {
593 cout << " 0\t";
594 }
595 }
596 // Geometry Collection is not a valid argument
597 catch(const IllegalArgumentException&) {
598 cout << " X\t";
599 }
600 catch(const std::exception& exc) {
601 cerr << exc.what() << endl;
602 }
603 }
604 cout << endl;
605 }
606
607 /////////////////////////////////////////////
608 // INTERSECTS
609 /////////////////////////////////////////////
610
611 cout << endl;
612 cout << " INTERSECTS ";
613 for(unsigned int i = 0; i < geoms->size(); i++) {
614 cout << "\t[" << i << "]";
615 }
616 cout << endl;
617 for(unsigned int i = 0; i < geoms->size(); i++) {
618 const Geometry* g1 = (*geoms)[i];
619 cout << " [" << i << "]\t";
620 for(unsigned int j = 0; j < geoms->size(); j++) {
621 const Geometry* g2 = (*geoms)[j];
622 try {
623 if(g1->intersects(g2)) {
624 cout << " 1\t";
625 }
626 else {
627 cout << " 0\t";
628 }
629 }
630 // Geometry Collection is not a valid argument
631 catch(const IllegalArgumentException&) {
632 cout << " X\t";
633 }
634 catch(const std::exception& exc) {
635 cerr << exc.what() << endl;
636 }
637 }
638 cout << endl;
639 }
640
641 /////////////////////////////////////////////
642 // CROSSES
643 /////////////////////////////////////////////
644
645 cout << endl;
646 cout << " CROSSES ";
647 for(unsigned int i = 0; i < geoms->size(); i++) {
648 cout << "\t[" << i << "]";
649 }
650 cout << endl;
651 for(unsigned int i = 0; i < geoms->size(); i++) {
652 const Geometry* g1 = (*geoms)[i];
653 cout << " [" << i << "]\t";
654 for(unsigned int j = 0; j < geoms->size(); j++) {
655 const Geometry* g2 = (*geoms)[j];
656 try {
657 if(g1->crosses(g2)) {
658 cout << " 1\t";
659 }
660 else {
661 cout << " 0\t";
662 }
663 }
664 // Geometry Collection is not a valid argument
665 catch(const IllegalArgumentException&) {
666 cout << " X\t";
667 }
668 catch(const std::exception& exc) {
669 cerr << exc.what() << endl;
670 }
671 }
672 cout << endl;
673 }
674
675 /////////////////////////////////////////////
676 // WITHIN
677 /////////////////////////////////////////////
678
679 cout << endl;
680 cout << " WITHIN ";
681 for(unsigned int i = 0; i < geoms->size(); i++) {
682 cout << "\t[" << i << "]";
683 }
684 cout << endl;
685 for(unsigned int i = 0; i < geoms->size(); i++) {
686 const Geometry* g1 = (*geoms)[i];
687 cout << " [" << i << "]\t";
688 for(unsigned int j = 0; j < geoms->size(); j++) {
689 const Geometry* g2 = (*geoms)[j];
690 try {
691 if(g1->within(g2)) {
692 cout << " 1\t";
693 }
694 else {
695 cout << " 0\t";
696 }
697 }
698 // Geometry Collection is not a valid argument
699 catch(const IllegalArgumentException&) {
700 cout << " X\t";
701 }
702 catch(const std::exception& exc) {
703 cerr << exc.what() << endl;
704 }
705 }
706 cout << endl;
707 }
708
709 /////////////////////////////////////////////
710 // CONTAINS
711 /////////////////////////////////////////////
712
713 cout << endl;
714 cout << " CONTAINS ";
715 for(unsigned int i = 0; i < geoms->size(); i++) {
716 cout << "\t[" << i << "]";
717 }
718 cout << endl;
719 for(unsigned int i = 0; i < geoms->size(); i++) {
720 const Geometry* g1 = (*geoms)[i];
721 cout << " [" << i << "]\t";
722 for(unsigned int j = 0; j < geoms->size(); j++) {
723 const Geometry* g2 = (*geoms)[j];
724 try {
725 if(g1->contains(g2)) {
726 cout << " 1\t";
727 }
728 else {
729 cout << " 0\t";
730 }
731 }
732 // Geometry Collection is not a valid argument
733 catch(const IllegalArgumentException&) {
734 cout << " X\t";
735 }
736 catch(const std::exception& exc) {
737 cerr << exc.what() << endl;
738 }
739 }
740 cout << endl;
741 }
742
743 /////////////////////////////////////////////
744 // OVERLAPS
745 /////////////////////////////////////////////
746
747 cout << endl;
748 cout << " OVERLAPS ";
749 for(unsigned int i = 0; i < geoms->size(); i++) {
750 cout << "\t[" << i << "]";
751 }
752 cout << endl;
753 for(unsigned int i = 0; i < geoms->size(); i++) {
754 const Geometry* g1 = (*geoms)[i];
755 cout << " [" << i << "]\t";
756 for(unsigned int j = 0; j < geoms->size(); j++) {
757 const Geometry* g2 = (*geoms)[j];
758 try {
759 if(g1->overlaps(g2)) {
760 cout << " 1\t";
761 }
762 else {
763 cout << " 0\t";
764 }
765 }
766 // Geometry Collection is not a valid argument
767 catch(const IllegalArgumentException&) {
768 cout << " X\t";
769 }
770 catch(const std::exception& exc) {
771 cerr << exc.what() << endl;
772 }
773 }
774 cout << endl;
775 }
776
777 /////////////////////////////////////////////
778 // RELATE
779 /////////////////////////////////////////////
780
781 cout << endl;
782 cout << " RELATE ";
783 for(unsigned int i = 0; i < geoms->size(); i++) {
784 cout << "\t[" << i << "]";
785 }
786 cout << endl;
787 for(unsigned int i = 0; i < geoms->size(); i++) {
788 const Geometry* g1 = (*geoms)[i];
789 cout << " [" << i << "]\t";
790 for(unsigned int j = 0; j < geoms->size(); j++) {
791 const Geometry* g2 = (*geoms)[j];
792 try {
793 // second argument is intersectionPattern
794 string pattern = "212101212";
795 if(g1->relate(g2, pattern)) {
796 cout << " 1\t";
797 }
798 else {
799 cout << " 0\t";
800 }
801
802 // get the intersectionMatrix itself
803 auto im = g1->relate(g2);
804 }
805 // Geometry Collection is not a valid argument
806 catch(const IllegalArgumentException&) {
807 cout << " X\t";
808 }
809 catch(const std::exception& exc) {
810 cerr << exc.what() << endl;
811 }
812 }
813 cout << endl;
814 }
815
816 /////////////////////////////////////////////
817 // EQUALS
818 /////////////////////////////////////////////
819
820 cout << endl;
821 cout << " EQUALS ";
822 for(unsigned int i = 0; i < geoms->size(); i++) {
823 cout << "\t[" << i << "]";
824 }
825 cout << endl;
826 for(unsigned int i = 0; i < geoms->size(); i++) {
827 const Geometry* g1 = (*geoms)[i];
828 cout << " [" << i << "]\t";
829 for(unsigned int j = 0; j < geoms->size(); j++) {
830 const Geometry* g2 = (*geoms)[j];
831 try {
832 if(g1->equals(g2)) {
833 cout << " 1\t";
834 }
835 else {
836 cout << " 0\t";
837 }
838 }
839 // Geometry Collection is not a valid argument
840 catch(const IllegalArgumentException&) {
841 cout << " X\t";
842 }
843 catch(const std::exception& exc) {
844 cerr << exc.what() << endl;
845 }
846 }
847 cout << endl;
848 }
849
850 /////////////////////////////////////////////
851 // EQUALS_EXACT
852 /////////////////////////////////////////////
853
854 cout << endl;
855 cout << "EQUALS_EXACT ";
856 for(unsigned int i = 0; i < geoms->size(); i++) {
857 cout << "\t[" << i << "]";
858 }
859 cout << endl;
860 for(unsigned int i = 0; i < geoms->size(); i++) {
861 const Geometry* g1 = (*geoms)[i];
862 cout << " [" << i << "]\t";
863 for(unsigned int j = 0; j < geoms->size(); j++) {
864 const Geometry* g2 = (*geoms)[j];
865 try {
866 // second argument is a tolerance
867 if(g1->equalsExact(g2, 0.5)) {
868 cout << " 1\t";
869 }
870 else {
871 cout << " 0\t";
872 }
873 }
874 // Geometry Collection is not a valid argument
875 catch(const IllegalArgumentException&) {
876 cout << " X\t";
877 }
878 catch(const std::exception& exc) {
879 cerr << exc.what() << endl;
880 }
881 }
882 cout << endl;
883 }
884
885 /////////////////////////////////////////////
886 // IS_WITHIN_DISTANCE
887 /////////////////////////////////////////////
888
889 cout << endl;
890 cout << "IS_WITHIN_DIST";
891 for(unsigned int i = 0; i < geoms->size(); i++) {
892 cout << "\t[" << i << "]";
893 }
894 cout << endl;
895 for(unsigned int i = 0; i < geoms->size(); i++) {
896 const Geometry* g1 = (*geoms)[i];
897 cout << " [" << i << "]\t";
898 for(unsigned int j = 0; j < geoms->size(); j++) {
899 const Geometry* g2 = (*geoms)[j];
900 try {
901 // second argument is the distance
902 if(g1->isWithinDistance(g2, 2)) {
903 cout << " 1\t";
904 }
905 else {
906 cout << " 0\t";
907 }
908 }
909 // Geometry Collection is not a valid argument
910 catch(const IllegalArgumentException&) {
911 cout << " X\t";
912 }
913 catch(const std::exception& exc) {
914 cerr << exc.what() << endl;
915 }
916 }
917 cout << endl;
918 }
919
920 #endif // RELATIONAL_OPERATORS
921
922 #if COMBINATIONS
923
924 ////////////////////////////////////////////////////////////////////////
925 // COMBINATIONS
926 ////////////////////////////////////////////////////////////////////////
927
928 cout << endl;
929 cout << "-------------------------------------------------------------------------------" << endl;
930 cout << "COMBINATIONS" << endl;
931 cout << "-------------------------------------------------------------------------------" << endl;
932
933 /////////////////////////////////////////////
934 // UNION
935 /////////////////////////////////////////////
936
937 // Make unions of all geoms
938 newgeoms = new vector<const Geometry*>;
939 for(unsigned int i = 0; i < geoms->size() - 1; i++) {
940 const Geometry* g1 = (*geoms)[i];
941 for(unsigned int j = i + 1; j < geoms->size(); j++) {
942 const Geometry* g2 = (*geoms)[j];
943 try {
944 Geometry* g3 = g1->Union(g2).release();
945 newgeoms->push_back(g3);
946 }
947 // It's illegal to union a collection ...
948 catch(const IllegalArgumentException&) {
949 //cerr <<ill.toString()<<"\n";
950 }
951 catch(const std::exception& exc) {
952 cerr << exc.what() << endl;
953 }
954 }
955 }
956
957 // Print all unions
958 cout << endl << "----- AND HERE ARE SOME UNION COMBINATIONS ------" << endl;
959 wkt_print_geoms(newgeoms);
960
961 // Delete the resulting geoms
962 for(unsigned int i = 0; i < newgeoms->size(); i++) {
963 delete(*newgeoms)[i];
964 }
965 delete newgeoms;
966
967
968 /////////////////////////////////////////////
969 // INTERSECTION
970 /////////////////////////////////////////////
971
972 // Compute intersection of adjacent geometries
973 newgeoms = new vector<const Geometry*>;
974 for(unsigned int i = 0; i < geoms->size() - 1; i++) {
975 const Geometry* g1 = (*geoms)[i];
976 for(unsigned int j = i + 1; j < geoms->size(); j++) {
977 const Geometry* g2 = (*geoms)[j];
978 try {
979 Geometry* g3 = g1->intersection(g2).release();
980 newgeoms->push_back(g3);
981 }
982 // Collection are illegal as intersection argument
983 catch(const IllegalArgumentException&) {
984 //cerr <<ill.toString()<<"\n";
985 }
986 catch(const std::exception& exc) {
987 cerr << exc.what() << endl;
988 }
989 }
990 }
991
992 cout << endl << "----- HERE ARE SOME INTERSECTIONS COMBINATIONS ------" << endl;
993 wkt_print_geoms(newgeoms);
994
995 // Delete the resulting geoms
996 for(unsigned int i = 0; i < newgeoms->size(); i++) {
997 delete(*newgeoms)[i];
998 }
999 delete newgeoms;
1000
1001 /////////////////////////////////////////////
1002 // DIFFERENCE
1003 /////////////////////////////////////////////
1004
1005 // Compute difference of adhiacent geometries
1006 newgeoms = new vector<const Geometry*>;
1007 for(unsigned int i = 0; i < geoms->size() - 1; i++) {
1008 const Geometry* g1 = (*geoms)[i];
1009 for(unsigned int j = i + 1; j < geoms->size(); j++) {
1010 const Geometry* g2 = (*geoms)[j];
1011 try {
1012 Geometry* g3 = g1->difference(g2).release();
1013 newgeoms->push_back(g3);
1014 }
1015 // Collection are illegal as difference argument
1016 catch(const IllegalArgumentException&) {
1017 //cerr <<ill.toString()<<"\n";
1018 }
1019 catch(const std::exception& exc) {
1020 cerr << exc.what() << endl;
1021 }
1022 }
1023 }
1024
1025 cout << endl << "----- HERE ARE SOME DIFFERENCE COMBINATIONS ------" << endl;
1026 wkt_print_geoms(newgeoms);
1027
1028 // Delete the resulting geoms
1029 for(unsigned int i = 0; i < newgeoms->size(); i++) {
1030 delete(*newgeoms)[i];
1031 }
1032 delete newgeoms;
1033
1034 /////////////////////////////////////////////
1035 // SYMMETRIC DIFFERENCE
1036 /////////////////////////////////////////////
1037
1038 // Compute symmetric difference of adhiacent geometries
1039 newgeoms = new vector<const Geometry*>;
1040 for(unsigned int i = 0; i < geoms->size() - 1; i++) {
1041 const Geometry* g1 = (*geoms)[i];
1042 for(unsigned int j = i + 1; j < geoms->size(); j++) {
1043 const Geometry* g2 = (*geoms)[j];
1044 try {
1045 Geometry* g3 = g1->symDifference(g2).release();
1046 newgeoms->push_back(g3);
1047 }
1048 // Collection are illegal as symdifference argument
1049 catch(const IllegalArgumentException&) {
1050 //cerr <<ill.toString()<<"\n";
1051 }
1052 catch(const std::exception& exc) {
1053 cerr << exc.what() << endl;
1054 }
1055 }
1056 }
1057
1058 cout << endl << "----- HERE ARE SYMMETRIC DIFFERENCES ------" << endl;
1059 wkt_print_geoms(newgeoms);
1060
1061 // Delete the resulting geoms
1062 for(unsigned int i = 0; i < newgeoms->size(); i++) {
1063 delete(*newgeoms)[i];
1064 }
1065 delete newgeoms;
1066
1067 #endif // COMBINATIONS
1068
1069 #if LINEMERGE
1070
1071 /////////////////////////////////////////////
1072 // LINEMERGE
1073 /////////////////////////////////////////////
1074 LineMerger lm;
1075 lm.add(geoms);
1076 auto mls = lm.getMergedLineStrings();
1077 newgeoms = new vector<const Geometry*>;
1078 for(unsigned int i = 0; i < mls.size(); i++) {
1079 newgeoms->push_back(mls[i].release());
1080 }
1081
1082 cout << endl << "----- HERE IS THE LINEMERGE OUTPUT ------" << endl;
1083 wkt_print_geoms(newgeoms);
1084
1085 // Delete the resulting geoms
1086 for(unsigned int i = 0; i < newgeoms->size(); i++) {
1087 delete(*newgeoms)[i];
1088 }
1089 delete newgeoms;
1090
1091 #endif // LINEMERGE
1092
1093 #if POLYGONIZE
1094
1095 /////////////////////////////////////////////
1096 // POLYGONIZE
1097 /////////////////////////////////////////////
1098 Polygonizer plgnzr;
1099 plgnzr.add(geoms);
1100 auto polys = plgnzr.getPolygons();
1101 newgeoms = new vector<const Geometry*>;
1102 for(unsigned int i = 0; i < polys.size(); i++) {
1103 newgeoms->push_back(polys[i].release());
1104 }
1105
1106 cout << endl << "----- HERE IS POLYGONIZE OUTPUT ------" << endl;
1107 wkt_print_geoms(newgeoms);
1108
1109 // Delete the resulting geoms
1110 for(unsigned int i = 0; i < newgeoms->size(); i++) {
1111 delete(*newgeoms)[i];
1112 }
1113 delete newgeoms;
1114
1115 #endif // POLYGONIZE
1116
1117 /////////////////////////////////////////////
1118 // CLEANUP
1119 /////////////////////////////////////////////
1120
1121 // Delete base geometries
1122 for(unsigned int i = 0; i < geoms->size(); i++) {
1123 delete(*geoms)[i];
1124 }
1125 delete geoms;
1126 }
1127
1128 int
main()1129 main()
1130 {
1131 cout << "GEOS " << geosversion() << " ported from JTS " << jtsport() << endl;
1132 try {
1133 do_all();
1134 }
1135 // All exception thrown by GEOS are subclasses of this
1136 // one, so this is a catch-all
1137 catch(const GEOSException& exc) {
1138 cerr << "GEOS Exception: " << exc.what() << "\n";
1139 exit(1);
1140 }
1141 catch(const exception& e) {
1142 cerr << "Standard exception thrown: " << e.what() << endl;
1143 exit(1);
1144 }
1145 // and this is a catch-all non standard ;)
1146 catch(...) {
1147 cerr << "unknown exception trown!\n";
1148 exit(1);
1149 }
1150
1151 // Unload is no more necessary
1152 //io::Unload::Release();
1153
1154 exit(0);
1155 }
1156