source: subversion/applications/utils/osm-extract/polygons/intersectpoly.cc @ 29220

Last change on this file since 29220 was 29220, checked in by frederik, 7 years ago

allow negative longitudes...

File size: 2.7 KB
Line 
1/*
2 * This little program computes the intersection between two .poly files,
3 * read on stdin, and writes a new .poly file on stdout.
4 *
5 * usage:
6 * cat p1.poly p2.poly ./intersectpoly > px.poly
7 *
8 * It can handle multipolygons but not holes; strange input data will
9 * yield unpredictable results (segfaults, most likely). Feel free to
10 * add appropriate error checking.
11 *
12 * compile with:
13 * c++ -o intersectpoly intersectpoly.cc -lgeos
14 *
15 * Written by Frederik Ramm <frederik@remote.org>, public domain.
16 */
17
18#include <stdio.h>
19#include <string.h>
20#include <stdlib.h>
21#include <vector>
22
23#include <geos/geom/LinearRing.h>
24#include <geos/geom/LinearRing.h>
25#include <geos/geom/Polygon.h>
26#include <geos/geom/MultiPolygon.h>
27#include <geos/geom/Coordinate.h>
28#include <geos/geom/CoordinateSequence.h>
29#include <geos/geom/CoordinateSequenceFactory.h>
30#include <geos/geom/GeometryFactory.h>
31
32using namespace geos::geom;
33using namespace std;
34
35GeometryFactory *gfactory;
36int count = 1;
37
38MultiPolygon* readpoly()
39{
40    int end = 0;
41    char buffer[256];
42    vector<Coordinate> *coords = new vector<Coordinate>();
43    vector<Geometry *> polys;
44
45    while (fgets(buffer, sizeof(buffer), stdin))
46    {
47        if (isspace(buffer[0]))
48        {
49            end = 0;
50            char *first = strpbrk(buffer, "0123456789.-");
51            char *second = strpbrk(first+1, " \t");
52            coords->push_back(Coordinate(atof(first), atof(second)));
53        }
54        else if (!strncmp(buffer, "END", 3))
55        {
56            if (end) break;
57            end++;
58            CoordinateSequence *cs = gfactory->getCoordinateSequenceFactory()->create(coords);
59            LinearRing *lr = gfactory->createLinearRing(cs);
60            Polygon *p = gfactory->createPolygon(lr, NULL);
61            polys.push_back(p);
62            coords = new vector<Coordinate>();
63        }
64        else
65        {
66            end = 0;
67        }
68    }
69    delete coords;
70    return gfactory->createMultiPolygon(polys);
71}
72
73void writepoly(const Geometry *p)
74{
75    printf("%d\n", count++);
76    CoordinateSequence *cs = p->getCoordinates();
77    for (int i=0; i<cs->getSize(); i++)
78    {
79        printf("   %10E   %10E\n", cs->getAt(i).x, cs->getAt(i).y);
80    }
81    printf("END\n");
82}
83
84int main(int argc, char **argv)
85{
86    gfactory = new GeometryFactory();
87    MultiPolygon *p1 = readpoly();
88    MultiPolygon *p2 = readpoly();
89    Geometry *intersect = p1->intersection(p2);
90    printf("none\n");
91    if (intersect->getGeometryTypeId() == GEOS_POLYGON)
92    {
93        writepoly(intersect);
94    }
95    else if (intersect->getGeometryTypeId() == GEOS_MULTIPOLYGON)
96    {
97        for (int i=0; i< intersect->getNumGeometries(); i++)
98        {
99            writepoly(intersect->getGeometryN(i));
100        }
101    }
102    printf("END\n");
103}
104
Note: See TracBrowser for help on using the repository browser.