/
osmium_noice.cpp
207 lines (167 loc) · 6.74 KB
/
osmium_noice.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
/*
osmium_noice
-----------------------------------------------------
extracts all non-icesheet polygons from an OSM file
for generating an antarctic icesheet polygon
by Christoph Hormann <chris_hormann@gmx.de>
based on osmium_toogr2 example
*/
#include <osmium/area/assembler.hpp>
#include <osmium/area/multipolygon_manager.hpp>
#include <osmium/geom/ogr.hpp>
#include <osmium/handler.hpp>
#include <osmium/handler/node_locations_for_ways.hpp>
#include <osmium/index/map/sparse_mem_array.hpp>
#include <osmium/index/node_locations_map.hpp>
#include <osmium/io/any_input.hpp>
#include <osmium/visitor.hpp>
#include <gdalcpp.hpp>
#include <getopt.h>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>
using index_type = osmium::index::map::SparseMemArray<osmium::unsigned_object_id_type, osmium::Location>;
using location_handler_type = osmium::handler::NodeLocationsForWays<index_type>;
template <class TProjection>
class MyOGRHandler : public osmium::handler::Handler {
gdalcpp::Layer m_layer_polygon;
osmium::geom::OGRFactory<TProjection>& m_factory;
std::size_t feature_count_noice = 0;
std::size_t feature_count_glacier = 0;
public:
MyOGRHandler(gdalcpp::Dataset& dataset, osmium::geom::OGRFactory<TProjection>& factory) :
m_layer_polygon(dataset, "noice", wkbMultiPolygon),
m_factory(factory) {
m_layer_polygon.add_field("id", OFTInteger, 10);
m_layer_polygon.add_field("type", OFTString, 30);
}
~MyOGRHandler() {
std::cout << "noice features converted: " << feature_count_noice << "\n";
std::cout << "glacier features converted: " << feature_count_glacier << "\n";
}
void area(const osmium::Area& area) {
const char* natural = area.tags()["natural"];
if (!natural) {
return;
}
// skip all areas that can exist on a glacier
// also ignore coastline tagging (since this is no area tag)
if (!std::strcmp(natural, "cliff") ||
!std::strcmp(natural, "sinkhole") ||
!std::strcmp(natural, "cave_entrance") ||
!std::strcmp(natural, "crevasse") ||
!std::strcmp(natural, "dune") ||
!std::strcmp(natural, "desert") ||
!std::strcmp(natural, "valley") ||
!std::strcmp(natural, "volcano") ||
!std::strcmp(natural, "coastline")) {
return;
}
const char* supraglacial = area.tags()["supraglacial"];
if (supraglacial && !std::strcmp(supraglacial, "yes")) {
return;
}
try {
gdalcpp::Feature feature(m_layer_polygon, m_factory.create_multipolygon(area));
feature.set_field("id", static_cast<double>(area.id()));
feature.set_field("type", natural);
feature.add_to_layer();
if (!std::strcmp(natural, "glacier")) {
++feature_count_glacier;
} else {
++feature_count_noice;
}
} catch (const osmium::geometry_error&) {
std::cerr << "Ignoring illegal geometry for area " << area.id() << " created from " << (area.from_way() ? "way" : "relation") << " with id=" << area.orig_id() << ".\n";
}
}
}; // class MyOGRHandler
/* ================================================== */
void print_help() {
std::cout << "osmium_noice [OPTIONS] [INFILE [OUTFILE]]\n\n"
<< "If INFILE is not given stdin is assumed.\n"
<< "If OUTFILE is not given 'ogr_out' is used.\n"
<< "\nOptions:\n"
<< " -h, --help This help message\n"
<< " -d, --debug Enable debug output\n"
<< " -f, --format=FORMAT Output OGR format (Default: 'SQLite')\n";
}
int main(int argc, char* argv[]) {
static struct option long_options[] = {
{"help", no_argument, nullptr, 'h'},
{"debug", no_argument, nullptr, 'd'},
{"format", required_argument, nullptr, 'f'},
{nullptr, 0, nullptr, 0}
};
std::string output_format{"SQLite"};
bool debug = false;
while (true) {
const int c = getopt_long(argc, argv, "hdf:", long_options, nullptr);
if (c == -1) {
break;
}
switch (c) {
case 'h':
print_help();
std::exit(0);
case 'd':
debug = true;
break;
case 'f':
output_format = optarg;
break;
default:
std::exit(1);
}
}
std::string input_filename;
std::string output_filename{"ogr_out"};
const int remaining_args = argc - optind;
if (remaining_args > 2) {
std::cerr << "Usage: " << argv[0] << " [OPTIONS] [INFILE [OUTFILE]]" << std::endl;
std::exit(1);
} else if (remaining_args == 2) {
input_filename = argv[optind];
output_filename = argv[optind+1];
} else if (remaining_args == 1) {
input_filename = argv[optind];
} else {
input_filename = "-";
}
osmium::io::File input_file{input_filename};
osmium::area::Assembler::config_type assembler_config;
if (debug) {
assembler_config.debug_level = 1;
}
osmium::area::MultipolygonManager<osmium::area::Assembler> mp_manager{assembler_config};
std::cerr << "Pass 1...\n";
osmium::relations::read_relations(input_file, mp_manager);
std::cerr << "Pass 1 done\n";
index_type index;
location_handler_type location_handler{index};
location_handler.ignore_errors();
osmium::geom::OGRFactory<> factory{};
CPLSetConfigOption("OGR_SQLITE_SYNCHRONOUS", "OFF");
gdalcpp::Dataset dataset{output_format, output_filename, gdalcpp::SRS{factory.proj_string()}, { "SPATIALITE=TRUE", "INIT_WITH_EPSG=no" }};
MyOGRHandler<decltype(factory)::projection_type> ogr_handler{dataset, factory};
std::cerr << "Pass 2...\n";
osmium::io::Reader reader{input_filename};
osmium::apply(reader, location_handler, ogr_handler, mp_manager.handler([&ogr_handler](const osmium::memory::Buffer& area_buffer) {
osmium::apply(area_buffer, ogr_handler);
}));
reader.close();
std::cerr << "Pass 2 done\n";
std::vector<osmium::object_id_type> incomplete_relations_ids;
mp_manager.for_each_incomplete_relation([&](const osmium::relations::RelationHandle& handle){
incomplete_relations_ids.push_back(handle->id());
});
if (!incomplete_relations_ids.empty()) {
std::cerr << "Warning! Some member ways missing for these multipolygon relations:";
for (const auto id : incomplete_relations_ids) {
std::cerr << " " << id;
}
std::cerr << "\n";
}
}