Osmium
0.1
|
00001 #ifndef OSMIUM_OSM_AREA_HPP 00002 #define OSMIUM_OSM_AREA_HPP 00003 00004 /* 00005 00006 Copyright 2011 Jochen Topf <jochen@topf.org> and others (see README). 00007 00008 This file is part of Osmium (https://github.com/joto/osmium). 00009 00010 Osmium is free software: you can redistribute it and/or modify it under the 00011 terms of the GNU Lesser General Public License or (at your option) the GNU 00012 General Public License as published by the Free Software Foundation, either 00013 version 3 of the Licenses, or (at your option) any later version. 00014 00015 Osmium is distributed in the hope that it will be useful, but WITHOUT ANY 00016 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 00017 PARTICULAR PURPOSE. See the GNU Lesser General Public License and the GNU 00018 General Public License for more details. 00019 00020 You should have received a copy of the Licenses along with Osmium. If not, see 00021 <http://www.gnu.org/licenses/>. 00022 00023 */ 00024 00025 #include <sys/types.h> 00026 #include <string.h> 00027 #include <stdio.h> 00028 #include <stdlib.h> 00029 #include <stdarg.h> 00030 #include <vector> 00031 #include <sstream> 00032 #include <iomanip> 00033 #include <map> 00034 00035 #ifdef OSMIUM_WITH_MULTIPOLYGON_PROFILING 00036 # include <osmium/utils/timer.h> 00037 # define START_TIMER(x) x##_timer.start(); 00038 # define STOP_TIMER(x) x##_timer.stop(); 00039 #else 00040 # define START_TIMER(x) 00041 # define STOP_TIMER(x) 00042 #endif // OSMIUM_WITH_MULTIPOLYGON_PROFILING 00043 00044 #ifdef OSMIUM_WITH_GEOS 00045 # include <geos/version.h> 00046 # include <geos/geom/Geometry.h> 00047 # include <geos/geom/Point.h> 00048 # include <geos/geom/LineString.h> 00049 # include <geos/geom/Polygon.h> 00050 # include <geos/geom/PrecisionModel.h> 00051 # include <geos/geom/CoordinateSequence.h> 00052 # include <geos/geom/CoordinateArraySequenceFactory.h> 00053 # include <geos/geom/LinearRing.h> 00054 # include <geos/geom/Polygon.h> 00055 # include <geos/geom/MultiPolygon.h> 00056 # include <geos/geom/MultiLineString.h> 00057 # include <geos/geom/prep/PreparedPolygon.h> 00058 # include <geos/io/WKTWriter.h> 00059 # include <geos/util/GEOSException.h> 00060 # include <geos/opLinemerge.h> 00061 # include <geos/operation/polygonize/Polygonizer.h> 00062 # include <geos/operation/distance/DistanceOp.h> 00063 # include <geos/opPolygonize.h> 00064 # include <geos/algorithm/LineIntersector.h> 00065 # include <geos/geomgraph/GeometryGraph.h> 00066 # include <geos/geomgraph/index/SegmentIntersector.h> 00067 00068 // this should come from /usr/include/geos/algorithm, but its missing there in some Ubuntu versions 00069 # include "../CGAlgorithms.h" 00070 #endif // OSMIUM_WITH_GEOS 00071 00072 #ifdef OSMIUM_WITH_SHPLIB 00073 # include <shapefil.h> 00074 #endif // OSMIUM_WITH_SHPLIB 00075 00076 #include <osmium/osm/way.hpp> 00077 #include <osmium/osm/relation.hpp> 00078 #include <osmium/geometry.hpp> 00079 00080 namespace Osmium { 00081 00082 namespace OSM { 00083 00084 enum innerouter_t { UNSET, INNER, OUTER }; 00085 enum direction_t { NO_DIRECTION, CLOCKWISE, COUNTERCLOCKWISE }; 00086 00087 #ifdef OSMIUM_WITH_GEOS 00088 class WayInfo { 00089 00090 friend class AreaFromRelation; 00091 00092 Osmium::OSM::Way *way; 00093 int used; 00094 int sequence; 00095 bool invert; 00096 bool duplicate; 00097 std::string errorhint; 00098 innerouter_t innerouter; 00099 innerouter_t orig_innerouter; 00100 geos::geom::Geometry *way_geom; 00101 int firstnode; 00102 int lastnode; 00103 bool tried; 00104 00105 WayInfo() { 00106 way = NULL; 00107 used = -1; 00108 innerouter = UNSET; 00109 orig_innerouter = UNSET; 00110 sequence = 0; 00111 invert = false; 00112 duplicate = false; 00113 way_geom = NULL; 00114 firstnode = -1; 00115 lastnode = -1; 00116 tried = false; 00117 } 00118 00119 WayInfo(Osmium::OSM::Way *w, innerouter_t io) { 00120 way = w; 00121 way_geom = w->create_geos_geometry(); 00122 orig_innerouter = io; 00123 used = -1; 00124 innerouter = UNSET; 00125 sequence = 0; 00126 invert = false; 00127 duplicate = false; 00128 tried = false; 00129 firstnode = w->get_first_node_id(); 00130 lastnode = w->get_last_node_id(); 00131 } 00132 00134 WayInfo(geos::geom::Geometry *g, int first, int last, innerouter_t io) { 00135 way = NULL; 00136 way_geom = g; 00137 orig_innerouter = io; 00138 used = -1; 00139 innerouter = UNSET; 00140 sequence = 0; 00141 invert = false; 00142 duplicate = false; 00143 tried = false; 00144 firstnode = first; 00145 lastnode = last; 00146 } 00147 00148 ~WayInfo() { 00149 delete way_geom; 00150 } 00151 00152 geos::geom::Point *get_firstnode_geom() { 00153 return (way ? way->get_first_node_geometry() : NULL); 00154 } 00155 00156 geos::geom::Point *get_lastnode_geom() { 00157 return (way ? way->get_last_node_geometry() : NULL); 00158 } 00159 00160 }; 00161 00162 class RingInfo { 00163 00164 friend class AreaFromRelation; 00165 00166 geos::geom::Polygon *polygon; 00167 direction_t direction; 00168 std::vector<WayInfo *> ways; 00169 std::vector<RingInfo *> inner_rings; 00170 bool nested; 00171 RingInfo *contained_by; 00172 int ring_id; 00173 00174 RingInfo() { 00175 polygon = NULL; 00176 direction = NO_DIRECTION; 00177 nested = false; 00178 contained_by = NULL; 00179 ring_id = -1; 00180 } 00181 }; 00182 #endif // OSMIUM_WITH_GEOS 00183 00185 class Area : public Object { 00186 00187 protected: 00188 00189 #ifdef OSMIUM_WITH_GEOS 00190 geos::geom::Geometry *geometry; 00191 00192 Area(geos::geom::Geometry *geom = NULL) : geometry(geom) { 00193 #else 00194 Area() { 00195 #endif // OSMIUM_WITH_GEOS 00196 # ifdef OSMIUM_WITH_JAVASCRIPT 00197 v8::HandleScope handle_scope; 00198 js_object_instance = v8::Persistent<v8::Object>::New(JavascriptTemplate::get<JavascriptTemplate>().create_instance(this)); 00199 # endif // OSMIUM_WITH_JAVASCRIPT 00200 } 00201 00202 ~Area() { 00203 #ifdef OSMIUM_WITH_GEOS 00204 delete(geometry); 00205 #endif // OSMIUM_WITH_GEOS 00206 } 00207 00208 public: 00209 00210 #ifdef OSMIUM_WITH_GEOS 00211 geos::geom::Geometry* get_geometry() const { 00212 return geometry; 00213 } 00214 #endif // OSMIUM_WITH_GEOS 00215 00216 #ifdef OSMIUM_WITH_JAVASCRIPT 00217 v8::Handle<v8::Value> js_from() const { 00218 const char *value = (get_type() == AREA_FROM_WAY) ? "way" : "relation"; 00219 return v8::String::NewSymbol(value); 00220 } 00221 00222 v8::Handle<v8::Value> js_geom() const; 00223 00224 struct JavascriptTemplate : public Osmium::OSM::Object::JavascriptTemplate { 00225 00226 JavascriptTemplate() : Osmium::OSM::Object::JavascriptTemplate() { 00227 js_template->SetAccessor(v8::String::NewSymbol("from"), accessor_getter<Area, &Area::js_from>); 00228 js_template->SetAccessor(v8::String::NewSymbol("geom"), accessor_getter<Area, &Area::js_geom>); 00229 } 00230 00231 }; 00232 00233 #endif // OSMIUM_WITH_JAVASCRIPT 00234 00235 }; // class Area 00236 00237 /*** 00238 * Area created from a way. 00239 * The way pointer given to the constructor will not be stored, all 00240 * needed attributes are copied. 00241 */ 00242 class AreaFromWay : public Area { 00243 00244 public: 00245 00246 WayNodeList m_node_list; 00247 00248 #ifdef OSMIUM_WITH_GEOS 00249 AreaFromWay(Way* way, geos::geom::Geometry* geom) : Area(geom) { 00250 #else 00251 AreaFromWay(Way* way) : Area() { 00252 #endif // OSMIUM_WITH_GEOS 00253 id(way->id()); 00254 version(way->version()); 00255 changeset(way->changeset()); 00256 timestamp(way->timestamp()); 00257 uid(way->uid()); 00258 user(way->user()); 00259 00260 tags(way->tags()); 00261 m_node_list = way->nodes(); 00262 } 00263 00264 ~AreaFromWay() { 00265 } 00266 00267 const WayNodeList& nodes() const { 00268 return m_node_list; 00269 } 00270 00271 WayNodeList& nodes() { 00272 return m_node_list; 00273 } 00274 00275 osm_object_type_t get_type() const { 00276 return AREA_FROM_WAY; 00277 } 00278 00279 }; // class AreaFromWay 00280 00281 /*** 00282 * Area created from a relation with tag type=multipolygon or type=boundary 00283 */ 00284 class AreaFromRelation : public Area { 00285 00286 bool boundary; 00287 00289 Relation *relation; 00290 00292 std::vector<Osmium::OSM::Way> member_ways; 00293 00295 int num_ways; 00296 00298 int missing_ways; 00299 00300 std::string geometry_error_message; 00301 00303 void (*callback)(Osmium::OSM::Area*); 00304 00306 bool attempt_repair; 00307 00308 bool ignore_tag(const std::string &s) { 00309 if (s=="type") return true; 00310 if (s=="created_by") return true; 00311 if (s=="source") return true; 00312 if (s=="note") return true; 00313 return false; 00314 } 00315 00316 bool same_tags(const Object* a, const Object* b) { 00317 if ((a == NULL) || (b == NULL)) return false; 00318 const TagList& at = a->tags(); 00319 const TagList& bt = b->tags(); 00320 std::map<std::string, std::string> aTags; 00321 00322 TagList::const_iterator end = at.end(); 00323 for (TagList::const_iterator it = at.begin(); it != end; ++it) { 00324 if (ignore_tag(it->key())) continue; 00325 aTags[it->key()] = it->value(); 00326 } 00327 end = bt.end(); 00328 for (TagList::const_iterator it = bt.begin(); it != end; ++it) { 00329 if (ignore_tag(it->key())) continue; 00330 if (aTags[it->key()] != it->value()) return false; 00331 aTags.erase(it->key()); 00332 } 00333 if (!aTags.empty()) return false; 00334 return true; 00335 } 00336 00338 bool merge_tags(Object* a, const Object* b) { 00339 bool rv = true; 00340 TagList& at = a->tags(); 00341 const TagList& bt = b->tags(); 00342 std::map<std::string, std::string> aTags; 00343 TagList::const_iterator end = at.end(); 00344 for (TagList::const_iterator it = at.begin(); it != end; ++it) { 00345 if (ignore_tag(it->key())) continue; 00346 aTags[it->key()] = it->value(); 00347 } 00348 end = bt.end(); 00349 for (TagList::const_iterator it = bt.begin(); it != end; ++it) { 00350 if (ignore_tag(it->key())) continue; 00351 if (aTags.find(it->key()) != aTags.end()) { 00352 if (aTags[it->key()] != it->value()) rv = false; 00353 } else { 00354 at.add(it->key(), it->value()); 00355 aTags[it->key()] = it->value(); 00356 } 00357 } 00358 return rv; 00359 } 00360 00361 bool untagged(const Object* r) { 00362 if (r == NULL) return true; 00363 const TagList& tags = r->tags(); 00364 if (tags.empty()) return true; 00365 TagList::const_iterator end = tags.end(); 00366 for (TagList::const_iterator it = tags.begin(); it != end; ++it) { 00367 if (! ignore_tag(it->key()) ) { 00368 return false; 00369 } 00370 } 00371 return true; 00372 } 00373 00374 public: 00375 00376 AreaFromRelation(Relation* r, bool b, int n, void (*callback)(Osmium::OSM::Area*), bool repair) : Area(), boundary(b), relation(r), callback(callback) { 00377 num_ways = n; 00378 missing_ways = n; 00379 #ifdef OSMIUM_WITH_GEOS 00380 geometry = NULL; 00381 #endif // OSMIUM_WITH_GEOS 00382 id(r->id()); 00383 attempt_repair = repair; 00384 } 00385 00386 #ifdef OSMIUM_WITH_MULTIPOLYGON_PROFILING 00387 static std::vector<std::pair<std::string, timer *> > timers; 00388 00389 static timer write_complex_poly_timer; 00390 static timer assemble_ways_timer; 00391 static timer assemble_nodes_timer; 00392 static timer make_one_ring_timer; 00393 static timer mor_polygonizer_timer; 00394 static timer mor_union_timer; 00395 static timer contains_timer; 00396 static timer extra_polygons_timer; 00397 static timer polygon_build_timer; 00398 static timer inner_ring_touch_timer; 00399 static timer polygon_intersection_timer; 00400 static timer multipolygon_build_timer; 00401 static timer multipolygon_write_timer; 00402 static timer error_write_timer; 00403 00404 static void init_timings() { 00405 timers.push_back(std::pair<std::string, timer *> (" thereof assemble_ways", &assemble_ways_timer)); 00406 timers.push_back(std::pair<std::string, timer *> (" thereof make_one_ring", &make_one_ring_timer)); 00407 timers.push_back(std::pair<std::string, timer *> (" thereof union", &mor_union_timer)); 00408 timers.push_back(std::pair<std::string, timer *> (" thereof polygonizer", &mor_polygonizer_timer)); 00409 timers.push_back(std::pair<std::string, timer *> (" thereof contains", &contains_timer)); 00410 timers.push_back(std::pair<std::string, timer *> (" thereof extra_polygons", &extra_polygons_timer)); 00411 timers.push_back(std::pair<std::string, timer *> (" thereof polygon_build", &polygon_build_timer)); 00412 timers.push_back(std::pair<std::string, timer *> (" thereof inner_ring_touch", &inner_ring_touch_timer)); 00413 timers.push_back(std::pair<std::string, timer *> (" thereof intersections", &polygon_intersection_timer)); 00414 timers.push_back(std::pair<std::string, timer *> (" thereof multipolygon_build", &multipolygon_build_timer)); 00415 timers.push_back(std::pair<std::string, timer *> (" thereof multipolygon_write", &multipolygon_write_timer)); 00416 timers.push_back(std::pair<std::string, timer *> (" thereof error_write", &error_write_timer)); 00417 } 00418 00419 static void print_timings() { 00420 for (unsigned int i=0; i<timers.size(); i++) { 00421 std::cerr << timers[i].first << ": " << *(timers[i].second) << std::endl; 00422 } 00423 } 00424 #endif // OSMIUM_WITH_MULTIPOLYGON_PROFILING 00425 00426 ~AreaFromRelation() { 00427 delete relation; 00428 member_ways.erase(member_ways.begin(), member_ways.end()); 00429 } 00430 00431 osm_object_type_t get_type() const { 00432 return AREA_FROM_RELATION; 00433 } 00434 00436 void add_member_way(Osmium::OSM::Way *way) { 00437 member_ways.push_back(*way); 00438 missing_ways--; 00439 } 00440 00442 bool is_complete() { 00443 return missing_ways == 0; 00444 } 00445 00446 void handle_complete_multipolygon() { 00447 #ifdef OSMIUM_WITH_GEOS 00448 if (! build_geometry()) { 00449 std::cerr << " geom build error: " << geometry_error_message << "\n"; 00450 } 00451 #endif // OSMIUM_WITH_GEOS 00452 callback(this); 00453 } 00454 00455 private: 00456 00468 #ifdef OSMIUM_WITH_GEOS 00469 geos::geom::LinearRing *create_non_intersecting_linear_ring(geos::geom::CoordinateSequence *orig_cs) { 00470 const std::vector<geos::geom::Coordinate>* coords = orig_cs->toVector(); 00471 int inv = coords->size(); 00472 int val = 0; 00473 int current = (inv + val) / 2; 00474 bool simple; 00475 00476 // find the longest non-intersecting stretch from the beginning 00477 // of the way. 00478 while (1) { 00479 std::vector<geos::geom::Coordinate> *vv = new std::vector<geos::geom::Coordinate>(coords->begin(), coords->begin() + current); 00480 geos::geom::CoordinateSequence *cs = geos::geom::CoordinateArraySequenceFactory::instance()->create(vv); 00481 geos::geom::LineString *a = Osmium::Geometry::geos_geometry_factory()->createLineString(cs); 00482 if (!(simple = a->isSimple())) { 00483 inv = current; 00484 } else { 00485 val = current; 00486 } 00487 delete a; 00488 if (current == (inv+val)/2) break; 00489 current = (inv + val) / 2; 00490 } 00491 00492 if (!simple) current--; 00493 00494 unsigned int cutoutstart = current; 00495 00496 inv = 0; 00497 val = coords->size(); 00498 current = (inv + val) / 2; 00499 00500 // find the longest non-intersecting stretch from the end 00501 // of the way. Note that this is likely to overlap with the 00502 // stretch found above - assume a 10-node way where nodes 3 00503 // and 7 are identical, then we will find the sequence 0..6 00504 // above, and 4..9 here! 00505 00506 while (1) { 00507 std::vector<geos::geom::Coordinate> *vv = new std::vector<geos::geom::Coordinate>(coords->begin() + current, coords->end()); 00508 geos::geom::CoordinateSequence *cs = geos::geom::CoordinateArraySequenceFactory::instance()->create(vv); 00509 geos::geom::LineString *a = Osmium::Geometry::geos_geometry_factory()->createLineString(cs); 00510 if (!(simple = a->isSimple())) { 00511 inv = current; 00512 } else { 00513 val = current; 00514 } 00515 delete a; 00516 if (current == (inv+val)/2) break; 00517 current = (inv + val) / 2; 00518 } 00519 if (!simple) current++; 00520 unsigned int cutoutend = current; 00521 00522 // assemble a new linear ring by cutting out the problematic bit. 00523 // if the "problematic bit" however is longer than half the way, 00524 // then try using the "problematic bit" by itself. 00525 00526 std::vector<geos::geom::Coordinate> *vv = new std::vector<geos::geom::Coordinate>(); 00527 if (cutoutstart < cutoutend) { 00528 unsigned int t = cutoutstart; 00529 cutoutstart = cutoutend; 00530 cutoutend = t; 00531 } 00532 if (cutoutstart-cutoutend > coords->size() / 2) { 00533 vv->insert(vv->end(), coords->begin() + cutoutend, coords->begin() + cutoutstart); 00534 vv->insert(vv->end(), vv->at(0)); 00535 } else { 00536 vv->insert(vv->end(), coords->begin(), coords->begin() + cutoutend); 00537 vv->insert(vv->end(), coords->begin() + cutoutstart, coords->end()); 00538 } 00539 geos::geom::CoordinateSequence *cs = geos::geom::CoordinateArraySequenceFactory::instance()->create(vv); 00540 geos::geom::LinearRing *a = Osmium::Geometry::geos_geometry_factory()->createLinearRing(cs); 00541 00542 // if this results in a valid ring, return it; else return NULL. 00543 00544 if (!a->isValid()) return NULL; 00545 geos::geom::LinearRing *b = dynamic_cast<geos::geom::LinearRing *>(a->clone()); 00546 //delete a; 00547 return b; 00548 } 00549 00557 RingInfo *make_one_ring(std::vector<WayInfo *> &ways, osm_object_id_t first, osm_object_id_t last, int ringcount, int sequence) { 00558 00559 // have we found a loop already? 00560 if (first && first == last) { 00561 geos::geom::CoordinateSequence *cs = geos::geom::CoordinateArraySequenceFactory::instance()->create((size_t)0, (size_t)0); 00562 geos::geom::LinearRing *lr = NULL; 00563 try { 00564 START_TIMER(mor_polygonizer); 00565 WayInfo **sorted_ways = new WayInfo*[sequence]; 00566 for (unsigned int i=0; i<ways.size(); i++) { 00567 if (ways[i]->used == ringcount) { 00568 sorted_ways[ways[i]->sequence] = ways[i]; 00569 } 00570 } 00571 for (int i=0; i<sequence; i++) { 00572 cs->add(dynamic_cast<geos::geom::LineString *>(sorted_ways[i]->way_geom)->getCoordinatesRO(), false, !sorted_ways[i]->invert); 00573 } 00574 delete[] sorted_ways; 00575 lr = Osmium::Geometry::geos_geometry_factory()->createLinearRing(cs); 00576 STOP_TIMER(mor_polygonizer); 00577 if (!lr->isSimple() || !lr->isValid()) { 00578 //delete lr; 00579 lr = NULL; 00580 if (attempt_repair) { 00581 lr = create_non_intersecting_linear_ring(cs); 00582 if (lr) { 00583 if (Osmium::debug()) 00584 std::cerr << "successfully repaired an invalid ring" << std::endl; 00585 } 00586 } 00587 if (!lr) return NULL; 00588 } 00589 bool ccw = geos::algorithm::CGAlgorithms::isCCW(lr->getCoordinatesRO()); 00590 RingInfo *rl = new RingInfo(); 00591 rl->direction = ccw ? COUNTERCLOCKWISE : CLOCKWISE; 00592 rl->polygon = Osmium::Geometry::geos_geometry_factory()->createPolygon(lr, NULL); 00593 return rl; 00594 } catch (const geos::util::GEOSException& exc) { 00595 if (Osmium::debug()) 00596 std::cerr << "Exception: " << exc.what() << std::endl; 00597 return NULL; 00598 } 00599 } 00600 00601 // have we not allocated anything yet, then simply start with first available way, 00602 // or return NULL if all are taken. 00603 if (!first) { 00604 for (unsigned int i=0; i<ways.size(); i++) { 00605 if (ways[i]->used != -1) continue; 00606 ways[i]->used = ringcount; 00607 ways[i]->sequence = 0; 00608 ways[i]->invert = false; 00609 RingInfo *rl = make_one_ring(ways, ways[i]->firstnode, ways[i]->lastnode, ringcount, 1); 00610 if (rl) { 00611 rl->ways.push_back(ways[i]); 00612 return rl; 00613 } 00614 ways[i]->used = -2; 00615 break; 00616 } 00617 return NULL; 00618 } 00619 00620 // try extending our current line at the rear end 00621 // since we are looking for a LOOP, no sense to try extending it at both ends 00622 // as we'll eventually get there anyway! 00623 00624 for (unsigned int i=0; i<ways.size(); i++) { 00625 if (ways[i]->used < 0) ways[i]->tried = false; 00626 } 00627 00628 for (unsigned int i=0; i<ways.size(); i++) { 00629 // ignore used ways 00630 if (ways[i]->used >= 0) continue; 00631 if (ways[i]->tried) continue; 00632 ways[i]->tried = true; 00633 00634 int old_used = ways[i]->used; 00635 if (ways[i]->firstnode == last) { 00636 // add way to end 00637 ways[i]->used = ringcount; 00638 ways[i]->sequence = sequence; 00639 ways[i]->invert = false; 00640 RingInfo *result = make_one_ring(ways, first, ways[i]->lastnode, ringcount, sequence+1); 00641 if (result) { 00642 result->ways.push_back(ways[i]); 00643 return result; 00644 } 00645 ways[i]->used = old_used; 00646 } else if (ways[i]->lastnode == last) { 00647 // add way to end, but turn it around 00648 ways[i]->used = ringcount; 00649 ways[i]->sequence = sequence; 00650 ways[i]->invert = true; 00651 RingInfo *result = make_one_ring(ways, first, ways[i]->firstnode, ringcount, sequence+1); 00652 if (result) { 00653 result->ways.push_back(ways[i]); 00654 return result; 00655 } 00656 ways[i]->used = old_used; 00657 } 00658 } 00659 // we have exhausted all combinations. 00660 return NULL; 00661 } 00662 00674 bool find_and_repair_holes_in_rings(std::vector<WayInfo *> *ways) { 00675 // collect the remaining debris (=unused ways) and find dangling nodes. 00676 00677 std::map<int,geos::geom::Point *> dangling_node_map; 00678 for (std::vector<WayInfo *>::iterator i(ways->begin()); i != ways->end(); i++) { 00679 if ((*i)->used < 0) { 00680 (*i)->innerouter = UNSET; 00681 (*i)->used = -1; 00682 for (int j=0; j<2; j++) { 00683 int nid = j ? (*i)->firstnode : (*i)->lastnode; 00684 if (dangling_node_map[nid]) { 00685 delete dangling_node_map[nid]; 00686 dangling_node_map[nid] = NULL; 00687 } else { 00688 dangling_node_map[nid] = j ? (*i)->get_firstnode_geom() : (*i)->get_lastnode_geom(); 00689 } 00690 } 00691 } 00692 } 00693 00694 do { 00695 int mindist_id = 0; 00696 double mindist = -1; 00697 int node1_id = 0; 00698 geos::geom::Point *node1 = NULL; 00699 geos::geom::Point *node2 = NULL; 00700 00701 // find one pair consisting of a random node from the list (node1) 00702 // plus the node that lies closest to it. 00703 for (std::map<int,geos::geom::Point *>::iterator i(dangling_node_map.begin()); i!= dangling_node_map.end(); i++) { 00704 if (!i->second) continue; 00705 if (node1 == NULL) { 00706 node1 = i->second; 00707 node1_id = i->first; 00708 i->second = NULL; 00709 mindist = -1; 00710 } else { 00711 # if GEOS_VERSION_MAJOR < 3 || (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR <= 2) 00712 double dist = geos::operation::distance::DistanceOp::distance(node1, (i->second)); // deprecated in newer version of GEOS 00713 # else 00714 double dist = geos::operation::distance::DistanceOp::distance(*node1, *(i->second)); 00715 # endif 00716 if ((dist < mindist) || (mindist < 0)) { 00717 mindist = dist; 00718 mindist_id = i->first; 00719 } 00720 } 00721 } 00722 00723 // if such a pair has been found, synthesize a connecting way. 00724 if (node1 && mindist > -1) { 00725 // if we find that there are dangling nodes but aren't 00726 // repairing - break out. 00727 if (!attempt_repair) return false; 00728 00729 // drop node2 from dangling map 00730 node2 = dangling_node_map[mindist_id]; 00731 dangling_node_map[mindist_id] = NULL; 00732 00733 std::vector<geos::geom::Coordinate> *c = new std::vector<geos::geom::Coordinate>; 00734 c->push_back(*(node1->getCoordinate())); 00735 c->push_back(*(node2->getCoordinate())); 00736 geos::geom::CoordinateSequence *cs = Osmium::Geometry::geos_geometry_factory()->getCoordinateSequenceFactory()->create(c); 00737 geos::geom::Geometry *geometry = (geos::geom::Geometry *) Osmium::Geometry::geos_geometry_factory()->createLineString(cs); 00738 ways->push_back(new WayInfo(geometry, node1_id, mindist_id, UNSET)); 00739 if (Osmium::debug()) 00740 std::cerr << "fill gap between nodes " << node1_id << " and " << mindist_id << std::endl; 00741 } else { 00742 break; 00743 } 00744 } while (1); 00745 00746 return true; 00747 } 00748 00749 00754 bool build_geometry() { 00755 std::vector<WayInfo *> ways; 00756 00757 // the timestamp of the multipolygon will be the maximum of the timestamp from the relation and from all member ways 00758 timestamp(relation->timestamp()); 00759 00760 // assemble all ways which are members of this relation into a 00761 // vector of WayInfo elements. this holds room for the way pointer 00762 // and some extra flags. 00763 00764 START_TIMER(assemble_ways); 00765 for (std::vector<Way>::iterator i(member_ways.begin()); i != member_ways.end(); i++) { 00766 if (i->timestamp() > timestamp()) { 00767 timestamp(i->timestamp()); 00768 } 00769 WayInfo *wi = new WayInfo(&(*i), UNSET); 00770 if (wi->way_geom) { 00771 geos::io::WKTWriter wkt; 00772 } else { 00773 delete wi; 00774 return geometry_error("invalid way geometry in multipolygon relation member"); 00775 } 00776 ways.push_back(wi); 00777 // TODO drop duplicate ways automatically in repair mode? 00778 // TODO maybe add INNER/OUTER instead of UNSET to enable later warnings on role mismatch 00779 } 00780 STOP_TIMER(assemble_ways); 00781 00782 std::vector<RingInfo *> ringlist; 00783 00784 // convenience defines to aid in clearing up on error return. 00785 #define clear_ringlist() \ 00786 for (std::vector<RingInfo *>::const_iterator rli(ringlist.begin()); rli != ringlist.end(); rli++) delete *rli; 00787 #define clear_wayinfo() \ 00788 for (std::vector<WayInfo *>::const_iterator win(ways.begin()); win != ways.end(); win++) delete *win; 00789 00790 // try and create as many closed rings as possible from the assortment 00791 // of ways. make_one_ring will automatically flag those that have been 00792 // used so they are not used again. 00793 00794 do { 00795 START_TIMER(make_one_ring); 00796 RingInfo *r = make_one_ring(ways, 0, 0, ringlist.size(), 0); 00797 STOP_TIMER(make_one_ring); 00798 if (r == NULL) break; 00799 r->ring_id = ringlist.size(); 00800 ringlist.push_back(r); 00801 } while (1); 00802 00803 if (ringlist.empty()) { 00804 // FIXME return geometry_error("no rings"); 00805 } 00806 00807 if (!find_and_repair_holes_in_rings(&ways)) { 00808 clear_ringlist(); 00809 clear_wayinfo(); 00810 return geometry_error("un-connectable dangling ends"); 00811 } 00812 00813 // re-run ring building, taking into account the newly created "repair" bits. 00814 // (in case there were no dangling bits, make_one_ring terminates quickly.) 00815 do { 00816 START_TIMER(make_one_ring); 00817 RingInfo *r = make_one_ring(ways, 0, 0, ringlist.size(), 0); 00818 STOP_TIMER(make_one_ring); 00819 if (r == NULL) break; 00820 r->ring_id = ringlist.size(); 00821 ringlist.push_back(r); 00822 } while (1); 00823 00824 if (ringlist.empty()) { 00825 clear_ringlist(); 00826 clear_wayinfo(); 00827 return geometry_error("no rings"); 00828 } 00829 00830 std::vector<geos::geom::Geometry *> *polygons = new std::vector<geos::geom::Geometry *>(); 00831 00832 geos::geom::MultiPolygon *mp = NULL; 00833 00834 // find out which ring contains which other ring, so we know 00835 // which are inner rings and which outer. don't trust the "role" 00836 // specifications. 00837 00838 START_TIMER(contains); 00839 00840 bool **contains = new bool*[ringlist.size()]; 00841 bool *contained_by_even_number = new bool[ringlist.size()]; 00842 00843 // reset array 00844 for (unsigned int i=0; i<ringlist.size(); i++) { 00845 contains[i] = new bool[ringlist.size()]; 00846 contained_by_even_number[i] = true; 00847 for (unsigned int j=0; j<ringlist.size(); j++) { 00848 contains[i][j] = false; 00849 } 00850 } 00851 00852 // build contains relationships. 00853 // we use contained_by_even_number as a helper for us to determine 00854 // whether something is an inner (false) or outer (true) ring. 00855 00856 for (unsigned int i=0; i<ringlist.size(); i++) { 00857 const geos::geom::prep::PreparedPolygon *pp = new geos::geom::prep::PreparedPolygon(ringlist[i]->polygon); 00858 for (unsigned int j=0; j<ringlist.size(); j++) { 00859 if (i==j) continue; 00860 if (contains[j][i]) continue; 00861 contains[i][j] = pp->contains(ringlist[j]->polygon); 00862 contained_by_even_number[j] ^= contains[i][j]; 00863 } 00864 delete pp; 00865 } 00866 00867 // we now have an array that has a true value whenever something is 00868 // contained by something else; if a contains b and b contains c, then 00869 // our array says that a contains b, b contains c, and a contains c. 00870 // thin out the array so that only direct relationships remain (and 00871 // the "a contains c" is dropped). 00872 00873 for (unsigned int i=0; i<ringlist.size(); i++) { 00874 for (unsigned j=0; j<ringlist.size(); j++) { 00875 if (contains[i][j]) { 00876 // see if there is an intermediary relationship 00877 for (unsigned int k=0; k<ringlist.size(); k++) { 00878 if (k==i) continue; 00879 if (k==j) continue; 00880 if (contains[i][k] && contains[k][j]) { 00881 // intermediary relationship exists; break this 00882 // one up. 00883 contains[i][j] = false; 00884 ringlist[j]->nested = true; 00885 break; 00886 } 00887 } 00888 } 00889 } 00890 } 00891 00892 // populate the "inner_rings" list and the "contained_by" pointer 00893 // in the ring list based on the data collected. the "contains" 00894 // array can be thrown away afterwards. 00895 00896 for (unsigned int i=0; i<ringlist.size(); i++) { 00897 for (unsigned int j=0; j<ringlist.size(); j++) { 00898 if (contains[i][j] && !contained_by_even_number[j]) { 00899 ringlist[j]->contained_by = ringlist[i]; 00900 ringlist[i]->inner_rings.push_back(ringlist[j]); 00901 } 00902 } 00903 delete[] contains[i]; 00904 } 00905 00906 delete[] contains; 00907 delete[] contained_by_even_number; 00908 STOP_TIMER(contains); 00909 00910 // now look at all enclosed (inner) rings that consist of only one way. 00911 // if such an inner ring has way tags, do the following: 00912 // * emit an extra polygon for the inner ring if the tags are different 00913 // from the relation's 00914 // * emit a warning, and ignore the inner ring, if the tags are the same 00915 // as for the relation 00916 00917 START_TIMER(extra_polygons); 00918 for (unsigned int i=0; i<ringlist.size(); i++) { 00919 if (ringlist[i]->contained_by) { 00920 if (ringlist[i]->ways.size() == 1 && !untagged(ringlist[i]->ways[0]->way)) { 00921 std::vector<geos::geom::Geometry *> *g = new std::vector<geos::geom::Geometry *>; 00922 if (ringlist[i]->direction == CLOCKWISE) { 00923 g->push_back(ringlist[i]->polygon->clone()); 00924 } else { 00925 geos::geom::LineString *tmp = dynamic_cast<geos::geom::LineString *>(ringlist[i]->polygon->getExteriorRing()->reverse()); 00926 geos::geom::LinearRing *reversed_ring = 00927 Osmium::Geometry::geos_geometry_factory()->createLinearRing(tmp->getCoordinates()); 00928 delete tmp; 00929 g->push_back(Osmium::Geometry::geos_geometry_factory()->createPolygon(reversed_ring, NULL)); 00930 } 00931 00932 geos::geom::MultiPolygon *special_mp = Osmium::Geometry::geos_geometry_factory()->createMultiPolygon(g); 00933 00934 if (same_tags(ringlist[i]->ways[0]->way, relation)) { 00935 // warning 00936 // warnings.insert("duplicate_tags_on_inner"); 00937 } else if (ringlist[i]->contained_by->ways.size() == 1 && same_tags(ringlist[i]->ways[0]->way, ringlist[i]->contained_by->ways[0]->way)) { 00938 // warning 00939 // warnings.insert("duplicate_tags_on_inner"); 00940 } else { 00941 Osmium::OSM::AreaFromWay *internal_mp = 00942 new Osmium::OSM::AreaFromWay(ringlist[i]->ways[0]->way, special_mp); 00943 callback(internal_mp); 00944 delete internal_mp; 00945 // AreaFromWay destructor deletes the 00946 // geometry, so avoid to delete it again. 00947 special_mp = NULL; 00948 } 00949 delete special_mp; 00950 } 00951 } 00952 } 00953 STOP_TIMER(extra_polygons); 00954 00955 // for all non-enclosed rings, assemble holes and build polygon. 00956 00957 START_TIMER(polygon_build) 00958 for (unsigned int i=0; i<ringlist.size(); i++) { 00959 // look only at outer, i.e. non-contained rings. each ends up as one polygon. 00960 if (ringlist[i] == NULL) continue; // can happen if ring has been deleted 00961 if (ringlist[i]->contained_by) continue; 00962 00963 std::vector<geos::geom::Geometry *> *holes = new std::vector<geos::geom::Geometry *>(); // ownership is later transferred to polygon 00964 00965 START_TIMER(inner_ring_touch) 00966 for (int j=0; j<((int)ringlist[i]->inner_rings.size()-1); j++) { 00967 if (!ringlist[i]->inner_rings[j]->polygon) continue; 00968 geos::geom::LinearRing *ring = (geos::geom::LinearRing *) ringlist[i]->inner_rings[j]->polygon->getExteriorRing(); 00969 00970 // check if some of the rings touch another ring. 00971 00972 for (unsigned int k=j + 1; k<ringlist[i]->inner_rings.size(); k++) { 00973 if (!ringlist[i]->inner_rings[k]->polygon) continue; 00974 const geos::geom::Geometry *compare = ringlist[i]->inner_rings[k]->polygon->getExteriorRing(); 00975 geos::geom::Geometry *inter = NULL; 00976 try { 00977 if (!ring->intersects(compare)) continue; 00978 inter = ring->intersection(compare); 00979 } catch (const geos::util::GEOSException& exc) { 00980 // nop; 00981 } 00982 if (inter && (inter->getGeometryTypeId() == geos::geom::GEOS_LINESTRING || inter->getGeometryTypeId() == geos::geom::GEOS_MULTILINESTRING)) { 00983 // touching inner rings 00984 // this is allowed, but we must fix them up into a valid 00985 // geometry 00986 geos::geom::Geometry *diff = ring->symDifference(compare); 00987 geos::operation::polygonize::Polygonizer *p = new geos::operation::polygonize::Polygonizer(); 00988 p->add(diff); 00989 std::vector<geos::geom::Polygon*>* polys = p->getPolygons(); 00990 if (polys && polys->size() == 1) { 00991 ringlist[i]->inner_rings[j]->polygon = polys->at(0); 00992 bool ccw = geos::algorithm::CGAlgorithms::isCCW(polys->at(0)->getExteriorRing()->getCoordinatesRO()); 00993 ringlist[i]->inner_rings[j]->direction = ccw ? COUNTERCLOCKWISE : CLOCKWISE; 00994 ringlist[i]->inner_rings[k]->polygon = NULL; 00995 j=-1; 00996 break; 00997 } 00998 } else { 00999 // other kind of intersect between inner rings; this is 01000 // not allwoed and will lead to an exception later when 01001 // building the MP 01002 } 01003 } 01004 } 01005 STOP_TIMER(inner_ring_touch) 01006 01007 for (unsigned int j=0; j<ringlist[i]->inner_rings.size(); j++) { 01008 if (!ringlist[i]->inner_rings[j]->polygon) continue; 01009 geos::geom::LinearRing *ring = (geos::geom::LinearRing *) ringlist[i]->inner_rings[j]->polygon->getExteriorRing(); 01010 01011 if (ringlist[i]->inner_rings[j]->direction == CLOCKWISE) { 01012 // reverse ring 01013 geos::geom::LineString *tmp = dynamic_cast<geos::geom::LineString *>(ring->reverse()); 01014 geos::geom::LinearRing *reversed_ring = 01015 Osmium::Geometry::geos_geometry_factory()->createLinearRing(tmp->getCoordinates()); 01016 delete tmp; 01017 holes->push_back(reversed_ring); 01018 } else { 01019 holes->push_back(ring); 01020 } 01021 } 01022 01023 geos::geom::LinearRing *ring = (geos::geom::LinearRing *) ringlist[i]->polygon->getExteriorRing(); 01024 if (ringlist[i]->direction == COUNTERCLOCKWISE) { 01025 geos::geom::LineString *tmp = dynamic_cast<geos::geom::LineString *>(ring->reverse()); 01026 geos::geom::LinearRing *reversed_ring = Osmium::Geometry::geos_geometry_factory()->createLinearRing(tmp->getCoordinates()); 01027 ring = reversed_ring; 01028 delete tmp; 01029 } else { 01030 ring = dynamic_cast<geos::geom::LinearRing *>(ring->clone()); 01031 } 01032 delete ringlist[i]->polygon; 01033 ringlist[i]->polygon = NULL; 01034 geos::geom::Polygon *p = NULL; 01035 bool valid = false; 01036 01037 try { 01038 p = Osmium::Geometry::geos_geometry_factory()->createPolygon(ring, holes); 01039 if (p) valid = p->isValid(); 01040 } catch (const geos::util::GEOSException& exc) { 01041 // nop 01042 if (Osmium::debug()) 01043 std::cerr << "Exception during creation of polygon for relation #" << relation->id() << ": " << exc.what() << " (treating as invalid polygon)" << std::endl; 01044 } 01045 if (!valid) { 01046 // polygon is invalid. 01047 clear_ringlist(); 01048 clear_wayinfo(); 01049 if (p) delete p; 01050 else delete ring; 01051 return geometry_error("invalid ring"); 01052 } else { 01053 polygons->push_back(p); 01054 for (unsigned int k=0; k<ringlist[i]->ways.size(); k++) { 01055 WayInfo *wi = ringlist[i]->ways[k]; 01056 // may have "hole filler" ways in there, not backed by 01057 // proper way and thus no tags: 01058 if (wi->way == NULL) continue; 01059 if (untagged(wi->way)) { 01060 // way not tagged - ok 01061 } else if (same_tags(relation, wi->way)) { 01062 // way tagged the same as relation/previous ways, ok 01063 } else if (untagged(relation)) { 01064 // relation untagged; use tags from way; ok 01065 merge_tags(relation, wi->way); 01066 } 01067 01068 wi->innerouter = OUTER; 01069 if (wi->orig_innerouter == INNER && wi->errorhint.empty()) { 01070 // warning: inner/outer mismatch 01071 } 01072 } 01073 // copy tags from relation into area 01074 tags(relation->tags()); 01075 } 01076 // later delete ringlist[i]; 01077 // ringlist[i] = NULL; 01078 } 01079 STOP_TIMER(polygon_build); 01080 01081 clear_ringlist(); 01082 clear_wayinfo(); 01083 if (polygons->empty()) { 01084 return geometry_error("no rings"); 01085 } 01086 01087 START_TIMER(multipolygon_build); 01088 bool valid = false; 01089 try { 01090 mp = Osmium::Geometry::geos_geometry_factory()->createMultiPolygon(polygons); 01091 valid = mp->isValid(); 01092 } catch (const geos::util::GEOSException& exc) { 01093 // nop 01094 }; 01095 STOP_TIMER(multipolygon_build); 01096 if (valid) { 01097 geometry = mp; 01098 return true; 01099 } 01100 return geometry_error("multipolygon invalid"); 01101 } 01102 01103 bool geometry_error(const char *message) { 01104 geometry_error_message = message; 01105 if (Osmium::debug()) { 01106 std::cerr << "building mp failed: " << geometry_error_message << std::endl; 01107 } 01108 geometry = NULL; 01109 return false; 01110 } 01111 01112 #endif // OSMIUM_WITH_GEOS 01113 01114 }; // class AreaFromRelation 01115 01116 } // namespace OSM 01117 01118 } // namespace Osmium 01119 01120 #endif // OSMIUM_OSM_AREA_HPP