[OSM-dev] libosmium and missing multipolygons

[Sylvain Melin]

Hi everyone.

I wrote a C++ program using libosmium to read an osm file, filter the 
data and extract the result to different target shapefiles.
I have a problem with some features not being written in shapefiles.

For example, I want to extract all water bodies in the bounding box 
latmin=44 latmax=45 lonmin=-1 lonmax=0.
In the resulting shapefile, the relation 2988 "Lit de la Dordogne" is 
missing : https://www.openstreetmap.org/relation/2988

The osmium::Area built from this relation seems to be empty.
area.num_rings() gives 0
area.tags() is empty

Trying to create a multipolygon from this area throws an 
osmium::geometry_error :
"Ignoring illegal geometry for area 5977 created from relation with id=2988.
area contains no rings (area_id=5977)"

Importing my osm file in QGis or converting it with ogr2ogr gives the 
appropriate result and the relation 2988 is correctly exported into the 
shapefile.

Do you have any idea of what's causing this problem and how it can be 
solved ?

Thank you for your help.

Here's the code of my osmium test program (it's a slight modification of 
https://github.com/osmcode/osm-gis-export/blob/master/src/osmium_toogr2.cpp) 
:

template <class TProjection>
class MyOGRHandler : public osmium::handler::Handler
{
     gdalcpp::Layer m_layer;

     osmium::geom::OGRFactory<TProjection>& m_factory;

public:

     MyOGRHandler(gdalcpp::Dataset& dataset, 
osmium::geom::OGRFactory<TProjection>& factory) :
         m_layer(dataset, "toto", wkbMultiPolygon),
         m_factory(factory)
     {
         m_layer.add_field("id", OFTReal, 10);
         m_layer.add_field("type", OFTString, 30);
     }

     void area(const osmium::Area& area)
     {
         try
         {
             gdalcpp::Feature feature(m_Output.m_layer, 
m_factory.create_multipolygon(area));
             feature.set_field("id", static_cast<double>(area.id()));
             feature.set_field("type", "toto");
             feature.add_to_layer();
         }
         catch (osmium::geometry_error& e)
         {
             std::cerr << "Ignoring illegal geometry for area "
                       << area.id()
                       << " created from "
                       << (area.from_way() ? "way" : "relation")
                       << " with id="
                       << area.orig_id() << ".\n";
             printf("%s\n\n",e.what());
         }
     }
};


int main(int argc, char* argv[])
{
     ...

             osmium::area::Assembler::config_type assembler_config;
             assembler_config.enable_debug_output(debug);
osmium::area::MultipolygonCollector<osmium::area::Assembler> 
collector(assembler_config);

             bool error = false;

             std::cerr << "Pass 1...\n";
             try
             {
                 osmium::io::Reader reader1(szIn);
                 collector.read_relations(reader1);
                 reader1.close();
                 printf("%s\n",collector.members_buffer().data());
             }
             catch(osmium::xml_error&)
             {
                 printf("Error reading file\n");
                 error = true;
             }
             std::cerr << "Pass 1 done\n";

             if(!error)
             {
                 index_type index_pos;
                 location_handler_type location_handler(index_pos);
                 location_handler.ignore_errors();

                 osmium::geom::OGRFactory<> factory {};

                 CPLSetConfigOption("OGR_SQLITE_SYNCHRONOUS", "FALSE");
                 gdalcpp::Dataset dataset{output_format, szOut, 
gdalcpp::SRS{factory.proj_string()}, { "SPATIALITE=TRUE" }};
MyOGRHandler<decltype(factory)::projection_type> ogr_handler(dataset, 
factory);

                 std::cerr << "Pass 2...\n";
                 osmium::io::Reader reader2(szIn);

                 osmium::apply(reader2, location_handler, ogr_handler, 
collector.handler([&ogr_handler](const osmium::memory::Buffer& 
area_buffer) {
                     osmium::apply(area_buffer, ogr_handler);
                 }));

                 reader2.close();
                 std::cerr << "Pass 2 done\n";

                 std::vector<const osmium::Relation*> 
incomplete_relations = collector.get_incomplete_relations();
                 if (!incomplete_relations.empty())
                 {
                     std::cerr << "Warning! Some member ways missing for 
these multipolygon relations:";
                     for (const auto* relation : incomplete_relations)
                     {
                         std::cerr << " " << relation->id();
                     }
                     std::cerr << "\n";
                 }
             }

}