source: subversion/applications/utils/export/osm2pgsql/build_geometry.cpp @ 15118

Last change on this file since 15118 was 14997, checked in by jonb, 11 years ago

Fix polygon ring directions using geos normalize()

File size: 11.2 KB
Line 
1/*
2#-----------------------------------------------------------------------------
3# Part of osm2pgsql utility
4#-----------------------------------------------------------------------------
5# By Artem Pavlenko, Copyright 2007
6#
7# This program is free software; you can redistribute it and/or
8# modify it under the terms of the GNU General Public License
9# as published by the Free Software Foundation; either version 2
10# of the License, or (at your option) any later version.
11#
12# This program is distributed in the hope that it will be useful,
13# but WITHOUT ANY WARRANTY; without even the implied warranty of
14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15# GNU General Public License for more details.
16#
17# You should have received a copy of the GNU General Public License
18# along with this program; if not, write to the Free Software
19# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
20#-----------------------------------------------------------------------------
21*/
22
23#include <iostream>
24#include <cstring>
25#include <cstdlib>
26
27/* Need to know which geos version we have to work out which headers to include */
28#include <geos/version.h>
29
30#if (GEOS_VERSION_MAJOR==3)
31/* geos trunk (3.0.0+) */
32#include <geos/geom/GeometryFactory.h>
33#include <geos/geom/CoordinateSequenceFactory.h>
34#include <geos/geom/Geometry.h>
35#include <geos/geom/LineString.h>
36#include <geos/geom/LinearRing.h>
37#include <geos/geom/MultiLineString.h>
38#include <geos/geom/Polygon.h>
39#include <geos/geom/Point.h>
40#include <geos/io/WKTReader.h>
41#include <geos/io/WKTWriter.h>
42#include <geos/util/GEOSException.h>
43#include <geos/opLinemerge.h>
44using namespace geos::geom;
45using namespace geos::io;
46using namespace geos::util;
47using namespace geos::operation::linemerge;
48#else
49/* geos-2.2.3 */
50#include <geos/geom.h>
51#include <geos/io.h>
52#include <geos/opLinemerge.h>
53using namespace geos;
54#endif
55
56#include "build_geometry.h"
57
58using namespace std;
59
60char *get_wkt_simple(osmNode *nodes, int count, int polygon, double *area, double *int_x, double *int_y) {
61    GeometryFactory gf;
62    auto_ptr<CoordinateSequence> coords(gf.getCoordinateSequenceFactory()->create(0, 2));
63
64    try
65    {
66        for (int i = 0; i < count ; i++) {
67            Coordinate c;
68            c.x = nodes[i].lon;
69            c.y = nodes[i].lat;
70            coords->add(c, 0);
71        }
72
73        auto_ptr<Geometry> geom;
74        if (polygon && (coords->getSize() >= 4) && (coords->getAt(coords->getSize() - 1).equals2D(coords->getAt(0)))) {
75            auto_ptr<LinearRing> shell(gf.createLinearRing(coords.release()));
76            geom = auto_ptr<Geometry>(gf.createPolygon(shell.release(),
77                                      new vector<Geometry *>));
78            geom->normalize(); // Fix direction of ring
79            *area = geom->getArea();
80            try {
81                std::auto_ptr<Point> pt(geom->getInteriorPoint());
82                *int_x = pt->getX();
83                *int_y = pt->getY();
84            } catch (...) {
85       // This happens on some unusual polygons, we'll ignore them for now
86       //std::cerr << std::endl << "Exception finding interior point" << std::endl;
87                *int_x = *int_y = 0.0;
88            }
89        } else {
90            *area = 0;
91            *int_x = *int_y = 0;
92            if (coords->getSize() < 2)
93                return NULL;
94            geom = auto_ptr<Geometry>(gf.createLineString(coords.release()));
95        }
96
97        WKTWriter wktw;
98        string wkt = wktw.write(geom.get());
99        return strdup(wkt.c_str());
100    }
101    catch (...)
102    {
103        std::cerr << std::endl << "excepton caught processing way" << std::endl;
104        return NULL;
105    }
106
107}
108
109
110typedef std::auto_ptr<Geometry> geom_ptr;
111
112struct Segment
113{
114      Segment(double x0_,double y0_,double x1_,double y1_)
115         :x0(x0_),y0(y0_),x1(x1_),y1(y1_) {}
116
117      double x0;
118      double y0;
119      double x1;
120      double y1;
121};
122
123struct Interior
124{
125        Interior(double x_, double y_)
126          : x(x_), y(y_) {}
127
128        Interior(Geometry *geom) {
129            try {
130                std::auto_ptr<Point> pt(geom->getInteriorPoint());
131                x = pt->getX();
132                y = pt->getY();
133            } catch (...) {
134                // This happens on some unusual polygons, we'll ignore them for now
135                //std::cerr << std::endl << "Exception finding interior point" << std::endl;
136                x=y=0.0;
137            }
138        }
139        double x, y;
140};
141
142static std::vector<std::string> wkts;
143static std::vector<Interior> interiors;
144static std::vector<double> areas;
145
146
147char * get_wkt(size_t index)
148{
149//   return wkts[index].c_str();
150        char *result;
151        result = (char*) std::malloc( wkts[index].length() + 1);
152        std::strcpy(result, wkts[index].c_str() );
153        return result;
154}
155
156void get_interior(size_t index, double *y, double *x)
157{
158        *x = interiors[index].x;
159        *y = interiors[index].y;
160}
161
162double get_area(size_t index)
163{
164    return areas[index];
165}
166
167void clear_wkts()
168{
169   wkts.clear();
170   interiors.clear();
171   areas.clear();
172}
173
174static int coords2nodes(CoordinateSequence * coords, struct osmNode ** nodes) {
175    size_t                      num_coords;
176    size_t                      i;
177    Coordinate          coord;
178
179    num_coords = coords->getSize();
180    *nodes = (struct osmNode *) malloc(num_coords * sizeof(struct osmNode));
181
182    for (i = 0; i < num_coords; i++) {
183        coord = coords->getAt(i);
184        (*nodes)[i].lon = coord.x;
185        (*nodes)[i].lat = coord.y;
186    }
187    return num_coords;
188}
189
190int parse_wkt(const char * wkt, struct osmNode *** xnodes, int ** xcount, int * polygon) {
191    GeometryFactory             gf;
192    WKTReader           reader(&gf);
193    std::string         wkt_string(wkt);
194    Geometry *          geometry;
195    const Geometry *    subgeometry;
196    GeometryCollection *        gc;
197    CoordinateSequence *        coords;
198    size_t                      num_geometries;
199    size_t                      i;
200       
201    *polygon = 0;
202    try {
203        geometry = reader.read(wkt_string);
204        switch (geometry->getGeometryTypeId()) {
205            // Single geometries
206            case GEOS_POLYGON:
207                // Drop through
208            case GEOS_LINEARRING:
209                *polygon = 1;
210                // Drop through
211            case GEOS_POINT:
212                // Drop through
213            case GEOS_LINESTRING:
214                *xnodes = (struct osmNode **) malloc(2 * sizeof(struct osmNode *));
215                *xcount = (int *) malloc(sizeof(int));
216                coords = geometry->getCoordinates();
217                (*xcount)[0] = coords2nodes(coords, &((*xnodes)[0]));
218                (*xnodes)[1] = NULL;
219                delete coords;
220                break;
221            // Geometry collections
222            case GEOS_MULTIPOLYGON:
223                *polygon = 1;
224                // Drop through
225            case GEOS_MULTIPOINT:
226                // Drop through
227            case GEOS_MULTILINESTRING:
228                gc = (GeometryCollection *) geometry;
229                num_geometries = gc->getNumGeometries();
230                *xnodes = (struct osmNode **) malloc((num_geometries + 1) * sizeof(struct osmNode *));
231                *xcount = (int *) malloc(num_geometries * sizeof(int));
232                for (i = 0; i < num_geometries; i++) {
233                    subgeometry = gc->getGeometryN(i);
234                    coords = subgeometry->getCoordinates();
235                    (*xcount)[0] = coords2nodes(coords, &((*xnodes)[i]));
236                    delete coords;
237                }
238                (*xnodes)[i] = NULL;
239                break;
240            default:
241                std::cerr << std::endl << "unexpected object type while processing PostGIS data" << std::endl;
242                delete geometry;
243                return -1;
244        }
245        delete geometry;
246    } catch (...) {
247        std::cerr << std::endl << "excepton caught parsing PostGIS data" << std::endl;
248        return -1;
249    }
250    return 0;
251}
252
253size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int make_polygon) {
254    size_t wkt_size = 0;
255    std::auto_ptr<std::vector<Geometry*> > lines(new std::vector<Geometry*>);
256    GeometryFactory gf;
257    auto_ptr<Geometry> geom;
258
259    try
260    {
261        for (int c=0; xnodes[c]; c++) {
262            auto_ptr<CoordinateSequence> coords(gf.getCoordinateSequenceFactory()->create(0, 2));
263            for (int i = 0; i < xcount[c]; i++) {
264                struct osmNode *nodes = xnodes[c];
265                Coordinate c;
266                c.x = nodes[i].lon;
267                c.y = nodes[i].lat;
268                coords->add(c, 0);
269            }
270            if (coords->getSize() > 1) {
271                geom = auto_ptr<Geometry>(gf.createLineString(coords.release()));
272                lines->push_back(geom.release());
273            }
274        }
275
276        //geom_ptr segment(0);
277        geom_ptr mline (gf.createMultiLineString(lines.release()));
278        //geom_ptr noded (segment->Union(mline.get()));
279        LineMerger merger;
280        //merger.add(noded.get());
281        merger.add(mline.get());
282        std::auto_ptr<std::vector<LineString *> > merged(merger.getMergedLineStrings());
283        WKTWriter writer;
284
285        std::auto_ptr<LinearRing> exterior;
286        std::auto_ptr<std::vector<Geometry*> > interior(new std::vector<Geometry*>);
287        double ext_area = 0.0;
288        for (unsigned i=0 ;i < merged->size(); ++i)
289        {
290            std::auto_ptr<LineString> pline ((*merged ) [i]);
291            if (make_polygon && pline->getNumPoints() > 3 && pline->isClosed())
292            {
293                std::auto_ptr<LinearRing> ring(gf.createLinearRing(pline->getCoordinates()));
294                std::auto_ptr<Polygon> poly(gf.createPolygon(gf.createLinearRing(pline->getCoordinates()),0));
295                double poly_area = poly->getArea();
296                if (poly_area > ext_area) {
297                    if (ext_area > 0.0)
298                        interior->push_back(exterior.release());
299                    ext_area = poly_area;
300                    exterior = ring;
301                        //std::cerr << "Found bigger ring, area(" << ext_area << ") " << writer.write(poly.get()) << std::endl;
302                } else {
303                    interior->push_back(ring.release());
304                        //std::cerr << "Found inner ring, area(" << poly->getArea() << ") " << writer.write(poly.get()) << std::endl;
305                }
306            } else {
307                        //std::cerr << "polygon(" << osm_id << ") is no good: points(" << pline->getNumPoints() << "), closed(" << pline->isClosed() << "). " << writer.write(pline.get()) << std::endl;
308                std::string text = writer.write(pline.get());
309                wkts.push_back(text);
310                interiors.push_back(Interior(0,0));
311                areas.push_back(0.0);
312                wkt_size++;
313            }
314        }
315
316        if (ext_area > 0.0) {
317            std::auto_ptr<Polygon> poly(gf.createPolygon(exterior.release(), interior.release()));
318            poly->normalize(); // Fix direction of rings
319            std::string text = writer.write(poly.get());
320                    //std::cerr << "Result: area(" << poly->getArea() << ") " << writer.write(poly.get()) << std::endl;
321            wkts.push_back(text);
322            interiors.push_back(Interior(poly.get()));
323            areas.push_back(ext_area);
324            wkt_size++;
325        }
326    }
327    catch (...)
328    {
329        std::cerr << std::endl << "excepton caught processing way id=" << osm_id << std::endl;
330        wkt_size = 0;
331    }
332    return wkt_size;
333}
Note: See TracBrowser for help on using the repository browser.