0> In article <51194252.7080...@web.de>,
0> WanMil <URL:mailto:wmgc...@web.de> ("Wanmil") wrote:
Wanmil> It might be possible. The problem to be solved is that the
Wanmil> garmin map geometry is not rectangular. So at the moment I
Wanmil> don't know how the area size of a polygon can be calculated in
Wanmil> a metric unit.
I did this by preprocessing in C++, which adds a "computed:area" tag to
each area (including multipolygons if TWOPASS is defined):
#include <osmium.hpp>
#include <osmium/storage/byid.hpp>
#include <osmium/handler/coordinates_for_ways.hpp>
#include <osmium/handler/multipolygon.hpp>
#include <osmium/geometry/multipolygon.hpp>
#include <geos/geom/Geometry.h>
#include <geos/geom/CoordinateFilter.h>
#include <geos/geom/Envelope.h>
#include <proj_api.h>
#include <unicode/coll.h>
#include <unicode/stsearch.h>
#include <cmath>
#include <vector>
#include <unordered_map>
#include <autosprintf.h>
static inline double deg2rad(double r) { return r * M_PI / 180; }
class ProjectFromWGS84 : public geos::geom::CoordinateFilter
{
projPJ source;
projPJ target;
public:
ProjectFromWGS84(const char *target_projection)
: geos::geom::CoordinateFilter(),
source(pj_init_plus("+proj=latlon +datum=WGS84")),
target(pj_init_plus(target_projection))
{
assert(source);
assert(target);
}
~ProjectFromWGS84()
{
pj_free(source);
pj_free(target);
}
virtual void filter_ro(const geos::geom::Coordinate *)
{
// not used
}
virtual void filter_rw(geos::geom::Coordinate *c) const
{
c->x *= DEG_TO_RAD;
c->y *= DEG_TO_RAD;
c->z *= DEG_TO_RAD;
pj_transform(source, target, 1, 1, &c->x, &c->y, &c->z);
}
};
typedef Osmium::Storage::SparseTable<Osmium::OSM::Position> storage_sparsetable_t;
typedef Osmium::Storage::Mmap<Osmium::OSM::Position> storage_mmap_t;
typedef Osmium::Handler::CoordinatesForWays<storage_sparsetable_t, storage_mmap_t> cfw_handler_t;
#if 1
class TagMap : public std::unordered_map<std::string,UnicodeString>
{
public:
TagMap(const Osmium::OSM::TagList& tags) {
reserve(tags.size());
for (const Osmium::OSM::Tag& tag: tags)
operator[](tag.key()) = UnicodeString::fromUTF8(tag.value());
}
operator Osmium::OSM::TagList() {
Osmium::OSM::TagList l;
for (const value_type& p: *this) {
std::string s;
p.second.toUTF8String(s);
l.add(p.first.c_str(), s.c_str());
}
return l;
}
bool contains(const key_type& k) const {
return find(k) != end();
}
const mapped_type& get(const key_type& k) const {
return contains(k) ? at(k) : empty_value;
}
private:
static const mapped_type empty_value;
};
const TagMap::mapped_type TagMap::empty_value = "";
template<typename Base>
class ModifyTagsHandler : public Base {
UErrorCode collator_error;
Collator *const collator;
public:
template<typename... Args>
ModifyTagsHandler(Args&&... args)
: Base(std::forward<Args>(args)...),
collator_error(U_ZERO_ERROR),
collator(Collator::createInstance(collator_error))
{
collator->setStrength(Collator::PRIMARY);
}
void node(const shared_ptr<Osmium::OSM::Node>& n) {
do_tag_update(n->tags(), n->get_type());
Base::node(n);
}
void way(const shared_ptr<Osmium::OSM::Way>& w) {
do_tag_update(w->tags(), w->get_type());
Base::way(w);
}
void relation(const shared_ptr<Osmium::OSM::Relation>& r) {
do_tag_update(r->tags(), r->get_type());
Base::relation(r);
}
private:
void do_tag_update(Osmium::OSM::TagList& tags, osm_object_type_t) {
using std::string;
TagMap m(tags);
UnicodeString oper = m.get("operator");
if (!oper.isEmpty()) {
static const Locale gb_locale("en", "GB");
UErrorCode status = U_ZERO_ERROR;
StringSearch iter(oper, m.get("name"), gb_locale, NULL, status);
//iter.getCollator()->setStrength(Collator::PRIMARY);
if (iter.first(status) != USEARCH_DONE) {
// name contains operator
m.erase("operator");
} else if (m.get("amenity") == "post_box" && oper.caseCompare("Royal Mail", U_FOLD_CASE_DEFAULT))
m.erase("operator");
else if (m.get("amenity") == "telephone" && oper.caseCompare("BT", U_FOLD_CASE_DEFAULT))
m.erase("operator");
}
tags = m;
}
};
#endif
#if 0
template<typename Callback>
class Multipolygon : public Base {
using std::vector;
/// a list of areas that need to be completed
vector<Osmium::OSM::AreaFromRelation*> m_areas;
// a map from way_id to a vector of indexes into the areas array
// this is used to find in which multipolygon relations a way is
typedef google::sparse_hash_map<osm_object_id_t, vector<osm_object_id_t> > way2areaidx_t;
way2areaidx_t m_way2areaidx;
bool m_attempt_repair;
Callback& m_callback_area;
uint64_t m_count_ways_in_all_areas;
public:
class Pass1Handler;
class Pass2Handler;
Multipolygon(bool attempt_repair, Callback& callback_area)
: Base(),
m_areas(),
m_way2areaidx(),
m_attempt_repair(attempt_repair),
m_callback_area(callback_area),
m_count_ways_in_all_areas(0) {
}
private:
class Pass1Handler {
void relation(const shared_ptr<Osmium::OSM::Relation const>& relation) {
const char* type = relation->tags().get_tag_by_key("type");
if (!type)
return;
const bool is_boundary = !strcmp(type, "boundary");
const bool is_multipolygon = !strcmp(type, "multipolygon");
if (!is_boundary && !is_multipolygon)
return;
int num_ways = 0;
for (const Osmium::OSM::RelationMember* member: relation->members()) {
if (member->type() == 'w') {
m_way2areaidx[member->ref()].push_back(m_areas.size());
++num_ways;
}
}
m_count_ways_in_all_areas += num_ways;
Osmium::OSM::AreaFromRelation* area = new Osmium::OSM::AreaFromRelation(new Osmium::OSM::Relation(*relation), is_boundary, num_ways, m_callback_area, m_attempt_repair);
m_areas.push_back(area);
}
};
class Pass2Handler {
void way(const shared_ptr<Osmium::OSM::Way>& way) {
way2areaidx_t::const_iterator way2areaidx_iterator(m_way2areaidx.find(way->id()));
if (way2areaidx_iterator == m_way2areaidx.end()) { // not in any relation
if (way->is_closed() && way->node_count() >= 4) { // way is closed and has enough nodes, build simple multipolygon
#ifdef OSMIUM_WITH_GEOS
Osmium::Geometry::Polygon polygon(*way);
Osmium::OSM::AreaFromWay* area = new Osmium::OSM::AreaFromWay(way.get(), polygon.create_geos_geometry());
#else
Osmium::OSM::AreaFromWay* area = new Osmium::OSM::AreaFromWay(way.get());
#endif // OSMIUM_WITH_GEOS
m_callback_area(area);
delete area;
}
return;
}
// is in at least one multipolygon relation
vector<osm_object_id_t>& v = way2areaidx_iterator->second;
for (const osm_object_id_t& i: v) {
Osmium::OSM::AreaFromRelation* area = m_areas[i];
if (!area) continue; // can't happen if pass 1 is correct
// store copy of current way in multipolygon
area->add_member_way(way.get());
if (area->is_complete()) {
area->handle_complete_multipolygon();
m_callback_area(area);
}
}
v.clear();
}
void after_ways() {
m_way2areaidx.clear();
}
};
};
#endif
template<typename Base>
class ComputeAreaHandler : public Base {
storage_sparsetable_t store_pos;
storage_mmap_t store_neg;
cfw_handler_t handler_cfw;
ProjectFromWGS84 project_north;
ProjectFromWGS84 project_world;
ProjectFromWGS84 project_south;
public:
template<typename... Args>
ComputeAreaHandler(Args&&... args)
: Base(std::forward<Args>(args)...),
handler_cfw(store_pos, store_neg),
project_north("+init=epsg:3408"),
project_world("+init=epsg:3410"),
project_south("+init=epsg:3409")
{
}
~ComputeAreaHandler() {
}
void init(Osmium::OSM::Meta& meta) {
handler_cfw.init(meta);
Base::init(meta);
}
void node(const shared_ptr<Osmium::OSM::Node>& n) {
handler_cfw.node(n);
Base::node(n);
}
void after_nodes() {
handler_cfw.after_nodes();
std::cerr << "nodes done" << std::endl;
Base::after_nodes();
}
void way(const shared_ptr<Osmium::OSM::Way>& w) {
#ifdef TWOPASS
handler_multipolygon->way(way);
#else
if (w->is_closed() && w->node_count() > 3) {
// remember, a triangle has 4 nodes (as one appears twice)
Osmium::Geometry::Polygon p(*w);
geos::geom::Geometry *g = p.create_geos_geometry();
addAreaTag(&*w, g);
delete g;
}
#endif
Base::way(w);
}
void tag_way(const shared_ptr<Osmium::OSM::Way>& w) {
if (w->is_closed() && w->node_count() > 3) {
// remember, a triangle has 4 nodes (as one appears twice)
Osmium::Geometry::Polygon p(*w);
geos::geom::Geometry *g = p.create_geos_geometry();
addAreaTag(&*w, g);
delete g;
}
}
void area(Osmium::OSM::Area *a)
{
geos::geom::Geometry *g = a->get_geometry();
addAreaTag(a, g);
Base::area(a);
}
private:
// Area of the supplied WGS-84 lon/lat geometry, in square metres
double computeArea(geos::geom::Geometry *g) {
if (!g)
return 0;
// We need to project to an equal-area projection for g->getArea() to give useful results
// We use azimuthal near the poles, and cylindrical near the equator, for precision
static const double max_lat = 60;
static const double hysteresis = 15;
static const geos::geom::Envelope north_rect(-180, 180, max_lat, 90);
static const geos::geom::Envelope world_rect(-180, 180, -max_lat, max_lat);
static const geos::geom::Envelope south_rect(-180, 180, -90, -60);
const geos::geom::Envelope *e = g->getEnvelopeInternal();
if (!e) return 0;
if (e->getMaxY() < (max_lat + hysteresis) && e->getMinY() > -(max_lat + hysteresis)) {
// most of the world => EPSG:3410
g->apply_rw(&project_world);
} else if (e->getMinY() > (max_lat - hysteresis)) {
// north polar region => EPSG:3408
g->apply_rw(&project_north);
} else if (e->getMaxY() < -(max_lat - hysteresis)) {
// south polar region => EPSG:3409
g->apply_rw(&project_south);
} else {
// big area - split into sub-areas
return computeArea(g->intersection(g->getFactory()->toGeometry(&north_rect)))
+ computeArea(g->intersection(g->getFactory()->toGeometry(&world_rect)))
+ computeArea(g->intersection(g->getFactory()->toGeometry(&south_rect)));
}
return g->getArea();
}
void addAreaTag(Osmium::OSM::Object *o, geos::geom::Geometry *g) {
std::ostringstream os;
os.precision(2);
os << std::fixed << computeArea(g);
o->tags().add("computed:area", os.str().c_str());
}
};
#if PJ_VERSION < 480
extern "C" {
void pj_clear_initcache(void);
}
#endif
int main(int argc, char *argv[])
{
std::string input;
std::string output;
if (argc > 3) {
std::cerr << "Usage: " << argv[0] << " [INFILE [OUTFILE]]" << std::endl;
exit(1);
}
if (argc > 1) {
input = argv[1];
}
if (argc > 2) {
output = argv[2];
}
Osmium::OSMFile outfile(output);
Osmium::Output::PBF out(outfile);
ModifyTagsHandler<ComputeAreaHandler<Osmium::Handler::Forward<decltype(out)>>> handler(&out);
//ComputeAreaHandler<Osmium::Handler::Forward<decltype(out)>> handler(&out);
//ComputeAreaHandler<Osmium::Handler::Forward<decltype(out)>> handler(&out);
Osmium::OSMFile infile(input);
infile.read(handler);
// clean up to help Valgrind
pj_clear_initcache();
xmlCleanupParser();
xmlCleanupCharEncodingHandlers();
return 0;
}
I'm stuck with an old osmium that can't handle 64-bit ids, but it should
still compile and operate with a more up-to-date version. I really
ought to get around to writing this up for the OSM wiki...
_______________________________________________
mkgmap-dev mailing list
mkgmap-dev@lists.mkgmap.org.uk
http://lists.mkgmap.org.uk/mailman/listinfo/mkgmap-dev
