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

Last change on this file since 12445 was 9441, checked in by andreas, 12 years ago

add missing include to compile with gcc-4.3

File size: 8.6 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/*
28  # the only reason we include this C header is to obtain the GEOS_VERSION_MAJOR
29  # definition to work out which headers we _really_ want. This version
30  # number is not present in any of the C++ headers (or at least it wasn't
31  # when I originally added those lines, probably 12+ months ago). This
32  # seemed the only easy way to support both geos versions.
33*/
34#include <geos_c.h>
35
36#if (GEOS_VERSION_MAJOR==3)
37/* geos trunk (3.0.0rc) */
38#include <geos/geom/GeometryFactory.h>
39#include <geos/geom/CoordinateSequenceFactory.h>
40#include <geos/geom/Geometry.h>
41#include <geos/geom/LineString.h>
42#include <geos/geom/LinearRing.h>
43#include <geos/geom/MultiLineString.h>
44#include <geos/geom/Polygon.h>
45#include <geos/geom/Point.h>
46#include <geos/io/WKTReader.h>
47#include <geos/io/WKTWriter.h>
48#include <geos/util/GEOSException.h>
49#include <geos/opLinemerge.h>
50using namespace geos::geom;
51using namespace geos::io;
52using namespace geos::util;
53using namespace geos::operation::linemerge;
54#else
55/* geos-2.2.3 */
56#include <geos/geom.h>
57#include <geos/io.h>
58#include <geos/opLinemerge.h>
59using namespace geos;
60#endif
61
62#include "build_geometry.h"
63
64using namespace std;
65
66char *get_wkt_simple(osmNode *nodes, int count, int polygon, double *area, double *int_x, double *int_y) {
67    GeometryFactory gf;
68    auto_ptr<CoordinateSequence> coords(gf.getCoordinateSequenceFactory()->create(0, 2));
69
70    try
71    {
72        for (int i = 0; i < count ; i++) {
73            Coordinate c;
74            c.x = nodes[i].lon;
75            c.y = nodes[i].lat;
76            coords->add(c, 0);
77        }
78
79        auto_ptr<Geometry> geom;
80        if (polygon && (coords->getSize() >= 4) && (coords->getAt(coords->getSize() - 1).equals2D(coords->getAt(0)))) {
81            auto_ptr<LinearRing> shell(gf.createLinearRing(coords.release()));
82            geom = auto_ptr<Geometry>(gf.createPolygon(shell.release(),
83                                      new vector<Geometry *>));
84            *area = geom->getArea();
85            try {
86                std::auto_ptr<Point> pt(geom->getInteriorPoint());
87                *int_x = pt->getX();
88                *int_y = pt->getY();
89            } catch (...) {
90       // This happens on some unusual polygons, we'll ignore them for now
91       //std::cerr << std::endl << "Exception finding interior point" << std::endl;
92                *int_x = *int_y = 0.0;
93            }
94        } else {
95            *area = 0;
96            *int_x = *int_y = 0;
97            if (coords->getSize() < 2)
98                return NULL;
99            geom = auto_ptr<Geometry>(gf.createLineString(coords.release()));
100        }
101
102        WKTWriter wktw;
103        string wkt = wktw.write(geom.get());
104        return strdup(wkt.c_str());
105    }
106    catch (...)
107    {
108        std::cerr << std::endl << "excepton caught processing way" << std::endl;
109        return NULL;
110    }
111
112}
113
114
115typedef std::auto_ptr<Geometry> geom_ptr;
116
117struct Segment
118{
119      Segment(double x0_,double y0_,double x1_,double y1_)
120         :x0(x0_),y0(y0_),x1(x1_),y1(y1_) {}
121
122      double x0;
123      double y0;
124      double x1;
125      double y1;
126};
127
128struct Interior
129{
130        Interior(double x_, double y_)
131          : x(x_), y(y_) {}
132
133        Interior(Geometry *geom) {
134            try {
135                std::auto_ptr<Point> pt(geom->getInteriorPoint());
136                x = pt->getX();
137                y = pt->getY();
138            } catch (...) {
139                // This happens on some unusual polygons, we'll ignore them for now
140                //std::cerr << std::endl << "Exception finding interior point" << std::endl;
141                x=y=0.0;
142            }
143        }
144        double x, y;
145};
146
147static std::vector<std::string> wkts;
148static std::vector<Interior> interiors;
149static std::vector<double> areas;
150
151
152char * get_wkt(size_t index)
153{
154//   return wkts[index].c_str();
155        char *result;
156        result = (char*) std::malloc( wkts[index].length() + 1);
157        std::strcpy(result, wkts[index].c_str() );
158        return result;
159}
160
161void get_interior(size_t index, double *y, double *x)
162{
163        *x = interiors[index].x;
164        *y = interiors[index].y;
165}
166
167double get_area(size_t index)
168{
169    return areas[index];
170}
171
172void clear_wkts()
173{
174   wkts.clear();
175   interiors.clear();
176   areas.clear();
177}
178
179size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int make_polygon) {
180    size_t wkt_size = 0;
181    std::auto_ptr<std::vector<Geometry*> > lines(new std::vector<Geometry*>);
182    GeometryFactory gf;
183    auto_ptr<Geometry> geom;
184
185    try
186    {
187        for (int c=0; xnodes[c]; c++) {
188            auto_ptr<CoordinateSequence> coords(gf.getCoordinateSequenceFactory()->create(0, 2));
189            for (int i = 0; i < xcount[c]; i++) {
190                struct osmNode *nodes = xnodes[c];
191                Coordinate c;
192                c.x = nodes[i].lon;
193                c.y = nodes[i].lat;
194                coords->add(c, 0);
195            }
196            if (coords->getSize() > 1) {
197                geom = auto_ptr<Geometry>(gf.createLineString(coords.release()));
198                lines->push_back(geom.release());
199            }
200        }
201
202        //geom_ptr segment(0);
203        geom_ptr mline (gf.createMultiLineString(lines.release()));
204        //geom_ptr noded (segment->Union(mline.get()));
205        LineMerger merger;
206        //merger.add(noded.get());
207        merger.add(mline.get());
208        std::auto_ptr<std::vector<LineString *> > merged(merger.getMergedLineStrings());
209        WKTWriter writer;
210
211        std::auto_ptr<LinearRing> exterior;
212        std::auto_ptr<std::vector<Geometry*> > interior(new std::vector<Geometry*>);
213        double ext_area = 0.0;
214        for (unsigned i=0 ;i < merged->size(); ++i)
215        {
216            std::auto_ptr<LineString> pline ((*merged ) [i]);
217            if (make_polygon && pline->getNumPoints() > 3 && pline->isClosed())
218            {
219                std::auto_ptr<LinearRing> ring(gf.createLinearRing(pline->getCoordinates()));
220                std::auto_ptr<Polygon> poly(gf.createPolygon(gf.createLinearRing(pline->getCoordinates()),0));
221                double poly_area = poly->getArea();
222                if (poly_area > ext_area) {
223                    if (ext_area > 0.0)
224                        interior->push_back(exterior.release());
225                    ext_area = poly_area;
226                    exterior = ring;
227                        //std::cerr << "Found bigger ring, area(" << ext_area << ") " << writer.write(poly.get()) << std::endl;
228                } else {
229                    interior->push_back(ring.release());
230                        //std::cerr << "Found inner ring, area(" << poly->getArea() << ") " << writer.write(poly.get()) << std::endl;
231                }
232            } else {
233                        //std::cerr << "polygon(" << osm_id << ") is no good: points(" << pline->getNumPoints() << "), closed(" << pline->isClosed() << "). " << writer.write(pline.get()) << std::endl;
234                std::string text = writer.write(pline.get());
235                wkts.push_back(text);
236                interiors.push_back(Interior(0,0));
237                areas.push_back(0.0);
238                wkt_size++;
239            }
240        }
241
242        if (ext_area > 0.0) {
243            std::auto_ptr<Polygon> poly(gf.createPolygon(exterior.release(), interior.release()));
244            std::string text = writer.write(poly.get());
245                    //std::cerr << "Result: area(" << poly->getArea() << ") " << writer.write(poly.get()) << std::endl;
246            wkts.push_back(text);
247            interiors.push_back(Interior(poly.get()));
248            areas.push_back(ext_area);
249            wkt_size++;
250        }
251    }
252    catch (...)
253    {
254        std::cerr << std::endl << "excepton caught processing way id=" << osm_id << std::endl;
255        wkt_size = 0;
256    }
257    return wkt_size;
258}
Note: See TracBrowser for help on using the repository browser.