source: subversion/applications/utils/export/osm2pgsql-intarray/build_geometry.cpp @ 28719

Last change on this file since 28719 was 25070, checked in by jonb, 9 years ago

osm2pgsql: Make use of prepared geometries from geos-3.1+. These dramatically speed up the polygon tests used by the advanced multipolygon handling code. Patch originally from Kompza.

File size: 18.7 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#include <exception>
27
28/* Need to know which geos version we have to work out which headers to include */
29#include <geos/version.h>
30
31/* geos (3.0.0+) */
32#if (GEOS_VERSION_MAJOR==3)
33#if (GEOS_VERSION_MINOR>=1)
34/* Prepared geometries are new in 3.1.0 */
35#define HAS_PREPARED_GEOMETRIES
36#include <geos/geom/prep/PreparedGeometryFactory.h>
37#include <geos/geom/prep/PreparedPolygon.h>
38#endif
39#include <geos/geom/GeometryFactory.h>
40#include <geos/geom/CoordinateSequenceFactory.h>
41#include <geos/geom/Geometry.h>
42#include <geos/geom/LineString.h>
43#include <geos/geom/LinearRing.h>
44#include <geos/geom/MultiLineString.h>
45#include <geos/geom/Polygon.h>
46#include <geos/geom/MultiPolygon.h>
47#include <geos/geom/Point.h>
48#include <geos/io/WKTReader.h>
49#include <geos/io/WKTWriter.h>
50#include <geos/util/GEOSException.h>
51#include <geos/opLinemerge.h>
52using namespace geos::geom;
53using namespace geos::io;
54using namespace geos::util;
55using namespace geos::operation::linemerge;
56#else
57/* geos-2.2.3 */
58#include <geos/geom.h>
59#include <geos/io.h>
60#include <geos/opLinemerge.h>
61using namespace geos;
62#endif
63
64#include "build_geometry.h"
65
66typedef std::auto_ptr<Geometry> geom_ptr;
67
68static std::vector<std::string> wkts;
69static std::vector<double> areas;
70
71
72char *get_wkt_simple(osmNode *nodes, int count, int polygon) {
73    GeometryFactory gf;
74    std::auto_ptr<CoordinateSequence> coords(gf.getCoordinateSequenceFactory()->create(0, 2));
75
76    try
77    {
78        for (int i = 0; i < count ; i++) {
79            Coordinate c;
80            c.x = nodes[i].lon;
81            c.y = nodes[i].lat;
82            coords->add(c, 0);
83        }
84
85        geom_ptr geom;
86        if (polygon && (coords->getSize() >= 4) && (coords->getAt(coords->getSize() - 1).equals2D(coords->getAt(0)))) {
87            std::auto_ptr<LinearRing> shell(gf.createLinearRing(coords.release()));
88            geom = geom_ptr(gf.createPolygon(shell.release(), new std::vector<Geometry *>));
89            geom->normalize(); // Fix direction of ring
90        } else {
91            if (coords->getSize() < 2)
92                return NULL;
93            geom = geom_ptr(gf.createLineString(coords.release()));
94        }
95
96        WKTWriter wktw;
97        std::string wkt = wktw.write(geom.get());
98        return strdup(wkt.c_str());
99    }
100    catch (std::bad_alloc)
101    {
102        std::cerr << std::endl << "Exception caught processing way. You are likelly running out of memory." << std::endl;
103        std::cerr << "Try in slim mode, using -s parameter." << std::endl;
104        return NULL;
105    }
106    catch (...)
107    {
108        std::cerr << std::endl << "Exception caught processing way" << std::endl;
109        return NULL;
110    }
111}
112
113
114size_t get_wkt_split(osmNode *nodes, int count, int polygon, double split_at) {
115    GeometryFactory gf;
116    std::auto_ptr<CoordinateSequence> coords(gf.getCoordinateSequenceFactory()->create(0, 2));
117    double area;
118    WKTWriter wktw;
119    size_t wkt_size = 0;
120
121    try
122    {
123        for (int i = 0; i < count ; i++) {
124            Coordinate c;
125            c.x = nodes[i].lon;
126            c.y = nodes[i].lat;
127            coords->add(c, 0);
128        }
129
130        geom_ptr geom;
131        if (polygon && (coords->getSize() >= 4) && (coords->getAt(coords->getSize() - 1).equals2D(coords->getAt(0)))) {
132            std::auto_ptr<LinearRing> shell(gf.createLinearRing(coords.release()));
133            geom = geom_ptr(gf.createPolygon(shell.release(), new std::vector<Geometry *>));
134            geom->normalize(); // Fix direction of ring
135            area = geom->getArea();
136            std::string wkt = wktw.write(geom.get());
137            wkts.push_back(wkt);
138            areas.push_back(area);
139            wkt_size++;
140        } else {
141            if (coords->getSize() < 2)
142                return 0;
143
144            double distance = 0;
145            std::auto_ptr<CoordinateSequence> segment;
146            segment = std::auto_ptr<CoordinateSequence>(gf.getCoordinateSequenceFactory()->create(0, 2));
147            segment->add(coords->getAt(0));
148            for(unsigned i=1; i<coords->getSize(); i++) {
149                segment->add(coords->getAt(i));
150                distance += coords->getAt(i).distance(coords->getAt(i-1));
151                if ((distance >= split_at) || (i == coords->getSize()-1)) {
152                    geom = geom_ptr(gf.createLineString(segment.release()));
153                    std::string wkt = wktw.write(geom.get());
154                    wkts.push_back(wkt);
155                    areas.push_back(0);
156                    wkt_size++;
157                    distance=0;
158                    segment = std::auto_ptr<CoordinateSequence>(gf.getCoordinateSequenceFactory()->create(0, 2));
159                    segment->add(coords->getAt(i));
160                }
161            }
162        }
163
164    }
165    catch (std::bad_alloc)
166    {
167        std::cerr << std::endl << "Exception caught processing way. You are likelly running out of memory." << std::endl;
168        std::cerr << "Try in slim mode, using -s parameter." << std::endl;
169        wkt_size = 0;
170    }
171    catch (...)
172    {
173        std::cerr << std::endl << "Exception caught processing way" << std::endl;
174        wkt_size = 0;
175    }
176    return wkt_size;
177}
178
179
180char * get_wkt(size_t index)
181{
182//   return wkts[index].c_str();
183        char *result;
184        result = (char*) std::malloc( wkts[index].length() + 1);
185        // At least give some idea of why we about to seg fault
186        if (!result) std::cerr << std::endl << "Unable to allocate memory: " << (wkts[index].length() + 1) << std::endl;
187        std::strcpy(result, wkts[index].c_str());
188        return result;
189}
190
191double get_area(size_t index)
192{
193    return areas[index];
194}
195
196void clear_wkts()
197{
198   wkts.clear();
199   areas.clear();
200}
201
202static int coords2nodes(CoordinateSequence * coords, struct osmNode ** nodes) {
203    size_t                      num_coords;
204    size_t                      i;
205    Coordinate          coord;
206
207    num_coords = coords->getSize();
208    *nodes = (struct osmNode *) malloc(num_coords * sizeof(struct osmNode));
209
210    for (i = 0; i < num_coords; i++) {
211        coord = coords->getAt(i);
212        (*nodes)[i].lon = coord.x;
213        (*nodes)[i].lat = coord.y;
214    }
215    return num_coords;
216}
217
218int parse_wkt(const char * wkt, struct osmNode *** xnodes, int ** xcount, int * polygon) {
219    GeometryFactory             gf;
220    WKTReader           reader(&gf);
221    std::string         wkt_string(wkt);
222    Geometry *          geometry;
223    const Geometry *    subgeometry;
224    GeometryCollection *        gc;
225    CoordinateSequence *        coords;
226    size_t                      num_geometries;
227    size_t                      i;
228       
229    *polygon = 0;
230    try {
231        geometry = reader.read(wkt_string);
232        switch (geometry->getGeometryTypeId()) {
233            // Single geometries
234            case GEOS_POLYGON:
235                // Drop through
236            case GEOS_LINEARRING:
237                *polygon = 1;
238                // Drop through
239            case GEOS_POINT:
240                // Drop through
241            case GEOS_LINESTRING:
242                *xnodes = (struct osmNode **) malloc(2 * sizeof(struct osmNode *));
243                *xcount = (int *) malloc(sizeof(int));
244                coords = geometry->getCoordinates();
245                (*xcount)[0] = coords2nodes(coords, &((*xnodes)[0]));
246                (*xnodes)[1] = NULL;
247                delete coords;
248                break;
249            // Geometry collections
250            case GEOS_MULTIPOLYGON:
251                *polygon = 1;
252                // Drop through
253            case GEOS_MULTIPOINT:
254                // Drop through
255            case GEOS_MULTILINESTRING:
256                gc = (GeometryCollection *) geometry;
257                num_geometries = gc->getNumGeometries();
258                *xnodes = (struct osmNode **) malloc((num_geometries + 1) * sizeof(struct osmNode *));
259                *xcount = (int *) malloc(num_geometries * sizeof(int));
260                for (i = 0; i < num_geometries; i++) {
261                    subgeometry = gc->getGeometryN(i);
262                    coords = subgeometry->getCoordinates();
263                    (*xcount)[i] = coords2nodes(coords, &((*xnodes)[i]));
264                    delete coords;
265                }
266                (*xnodes)[i] = NULL;
267                break;
268            default:
269                std::cerr << std::endl << "unexpected object type while processing PostGIS data" << std::endl;
270                delete geometry;
271                return -1;
272        }
273        delete geometry;
274    } catch (...) {
275        std::cerr << std::endl << "Exception caught parsing PostGIS data" << std::endl;
276        return -1;
277    }
278    return 0;
279}
280
281struct polygondata
282{
283    Polygon*        polygon;
284    LinearRing*     ring;
285    double          area;
286    int             iscontained;
287    unsigned        containedbyid;
288};
289
290static int polygondata_comparearea(const void* vp1, const void* vp2)
291{
292    const polygondata* p1 = (const polygondata*)vp1;
293    const polygondata* p2 = (const polygondata*)vp2;
294
295    if (p1->area == p2->area) return 0;
296    if (p1->area > p2->area) return -1;
297    return 1;
298}
299
300size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int make_polygon, int enable_multi, double split_at) {
301    size_t wkt_size = 0;
302    std::auto_ptr<std::vector<Geometry*> > lines(new std::vector<Geometry*>);
303    GeometryFactory gf;
304    geom_ptr geom;
305#ifdef HAS_PREPARED_GEOMETRIES
306    geos::geom::prep::PreparedGeometryFactory pgf;
307#endif
308
309    try
310    {
311        for (int c=0; xnodes[c]; c++) {
312            std::auto_ptr<CoordinateSequence> coords(gf.getCoordinateSequenceFactory()->create(0, 2));
313            for (int i = 0; i < xcount[c]; i++) {
314                struct osmNode *nodes = xnodes[c];
315                Coordinate c;
316                c.x = nodes[i].lon;
317                c.y = nodes[i].lat;
318                coords->add(c, 0);
319            }
320            if (coords->getSize() > 1) {
321                geom = geom_ptr(gf.createLineString(coords.release()));
322                lines->push_back(geom.release());
323            }
324        }
325
326        //geom_ptr segment(0);
327        geom_ptr mline (gf.createMultiLineString(lines.release()));
328        //geom_ptr noded (segment->Union(mline.get()));
329        LineMerger merger;
330        //merger.add(noded.get());
331        merger.add(mline.get());
332        std::auto_ptr<std::vector<LineString *> > merged(merger.getMergedLineStrings());
333        WKTWriter writer;
334
335        // Procces ways into lines or simple polygon list
336        polygondata* polys = new polygondata[merged->size()];
337
338        unsigned totalpolys = 0;
339        for (unsigned i=0 ;i < merged->size(); ++i)
340        {
341            std::auto_ptr<LineString> pline ((*merged ) [i]);
342            if (make_polygon && pline->getNumPoints() > 3 && pline->isClosed())
343            {
344                polys[totalpolys].polygon = gf.createPolygon(gf.createLinearRing(pline->getCoordinates()),0);
345                polys[totalpolys].ring = gf.createLinearRing(pline->getCoordinates());
346                polys[totalpolys].area = polys[totalpolys].polygon->getArea();
347                polys[totalpolys].iscontained = 0;
348                polys[totalpolys].containedbyid = 0;
349                if (polys[totalpolys].area > 0.0)
350                    totalpolys++;
351                else {
352                    delete(polys[totalpolys].polygon);
353                    delete(polys[totalpolys].ring);
354                }
355            }
356            else
357            {
358                        //std::cerr << "polygon(" << osm_id << ") is no good: points(" << pline->getNumPoints() << "), closed(" << pline->isClosed() << "). " << writer.write(pline.get()) << std::endl;
359                double distance = 0;
360                std::auto_ptr<CoordinateSequence> segment;
361                segment = std::auto_ptr<CoordinateSequence>(gf.getCoordinateSequenceFactory()->create(0, 2));
362                segment->add(pline->getCoordinateN(0));
363                for(unsigned i=1; i<pline->getNumPoints(); i++) {
364                    segment->add(pline->getCoordinateN(i));
365                    distance += pline->getCoordinateN(i).distance(pline->getCoordinateN(i-1));
366                    if ((distance >= split_at) || (i == pline->getNumPoints()-1)) {
367                        geom = geom_ptr(gf.createLineString(segment.release()));
368                        std::string wkt = writer.write(geom.get());
369                        wkts.push_back(wkt);
370                        areas.push_back(0);
371                        wkt_size++;
372                        distance=0;
373                        segment = std::auto_ptr<CoordinateSequence>(gf.getCoordinateSequenceFactory()->create(0, 2));
374                        segment->add(pline->getCoordinateN(i));
375                    }
376                }
377                //std::string text = writer.write(pline.get());
378                //wkts.push_back(text);
379                //areas.push_back(0.0);
380                //wkt_size++;
381            }
382        }
383
384        if (totalpolys)
385        {
386            qsort(polys, totalpolys, sizeof(polygondata), polygondata_comparearea);
387
388            unsigned toplevelpolygons = 0;
389            int istoplevelafterall;
390
391            for (unsigned i=0 ;i < totalpolys; ++i)
392            {
393                if (polys[i].iscontained != 0) continue;
394                toplevelpolygons++;
395#ifdef HAS_PREPARED_GEOMETRIES
396                const geos::geom::prep::PreparedGeometry* preparedtoplevelpolygon = pgf.create(polys[i].polygon);
397#endif
398
399                for (unsigned j=i+1; j < totalpolys; ++j)
400                {
401#ifdef HAS_PREPARED_GEOMETRIES
402                    // Does preparedtoplevelpolygon contain the smaller polygon[j]?
403                    if (polys[j].containedbyid == 0 && preparedtoplevelpolygon->contains(polys[j].polygon))
404#else
405                    // Does polygon[i] contain the smaller polygon[j]?
406                    if (polys[j].containedbyid == 0 && polys[i].polygon->contains(polys[j].polygon))
407#endif
408                    {
409                        // are we in a [i] contains [k] contains [j] situation
410                        // which would actually make j top level
411                        istoplevelafterall = 0;
412                        for (unsigned k=i+1; k < j; ++k)
413                        {
414                            if (polys[k].iscontained && polys[k].containedbyid == i && polys[k].polygon->contains(polys[j].polygon))
415                            {
416                                istoplevelafterall = 1;
417                                break;
418                            }
419#if 0
420                            else if (polys[k].polygon->intersects(polys[j].polygon) || polys[k].polygon->touches(polys[j].polygon))
421                            {
422                                // FIXME: This code does not work as intended
423                                // It should be setting the polys[k].ring in order to update this object
424                                // but the value of polys[k].polygon calculated is normally NULL
425
426                                // Add polygon this polygon (j) to k since they intersect
427                                // Mark ourselfs to be dropped (2), delete the original k
428                                Geometry* polyunion = polys[k].polygon->Union(polys[j].polygon);
429                                delete(polys[k].polygon);
430                                polys[k].polygon = dynamic_cast<Polygon*>(polyunion);
431                                polys[j].iscontained = 2; // Drop
432                                istoplevelafterall = 2;
433                                break;
434
435                            }
436#endif
437                        }
438                        if (istoplevelafterall == 0)
439                        {
440                            polys[j].iscontained = 1;
441                            polys[j].containedbyid = i;
442                        }
443                    }
444                }
445#ifdef HAS_PREPARED_GEOMETRIES
446                pgf.destroy(preparedtoplevelpolygon);
447#endif
448            }
449            // polys now is a list of ploygons tagged with which ones are inside each other
450
451            // List of polygons for multipolygon
452            std::auto_ptr<std::vector<Geometry*> > polygons(new std::vector<Geometry*>);
453
454            // For each top level polygon create a new polygon including any holes
455            for (unsigned i=0 ;i < totalpolys; ++i)
456            {
457                if (polys[i].iscontained != 0) continue;
458
459                // List of holes for this top level polygon
460                std::auto_ptr<std::vector<Geometry*> > interior(new std::vector<Geometry*>);
461                for (unsigned j=i+1; j < totalpolys; ++j)
462                {
463                   if (polys[j].iscontained == 1 && polys[j].containedbyid == i)
464                   {
465                       interior->push_back(polys[j].ring);
466                   }
467                }
468               
469                Polygon* poly(gf.createPolygon(polys[i].ring, interior.release()));
470                poly->normalize();
471                polygons->push_back(poly);
472            }
473
474            // Make a multipolygon if required
475            if ((toplevelpolygons > 1) && enable_multi)
476            {
477                std::auto_ptr<MultiPolygon> multipoly(gf.createMultiPolygon(polygons.release()));
478                std::string text = writer.write(multipoly.get());
479                wkts.push_back(text);
480                areas.push_back(multipoly->getArea());
481                wkt_size++;
482            }
483            else
484            {
485                for(unsigned i=0; i<toplevelpolygons; i++) 
486                {
487                    Polygon* poly = (Polygon*)polygons->at(i);
488                    std::string text = writer.write(poly);
489                    wkts.push_back(text);
490                    areas.push_back(poly->getArea());
491                    wkt_size++;
492                    delete(poly);
493                }
494            }
495        }
496
497        for (unsigned i=0; i < totalpolys; ++i)
498        {
499            delete(polys[i].polygon);
500        }
501        delete[](polys);
502    }
503    catch (std::exception& e)
504      {
505        std::cerr << std::endl << "Standard exception processing way_id "<< osm_id << ": " << e.what()  << std::endl;
506        wkt_size = 0;
507      }
508    catch (...)
509      {
510        std::cerr << std::endl << "Exception caught processing way id=" << osm_id << std::endl;
511        wkt_size = 0;
512      }
513
514    return wkt_size;
515}
Note: See TracBrowser for help on using the repository browser.