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";
}
}
}
_______________________________________________
dev mailing list
[email protected]
https://lists.openstreetmap.org/listinfo/dev
