[mkgmap-dev] Computing area (was: length () in relations file)
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::SparseTableOsmium::OSM::Position storage_sparsetable_t;
typedef Osmium::Storage::MmapOsmium::OSM::Position storage_mmap_t;
typedef Osmium::Handler::CoordinatesForWaysstorage_sparsetable_t, storage_mmap_t cfw_handler_t;
#if 1
class TagMap : public std::unordered_mapstd::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 = ;
templatetypename Base
class ModifyTagsHandler : public Base {
UErrorCode collator_error;
Collator *const collator;
public:
templatetypename... Args
ModifyTagsHandler(Args... args)
: Base(std::forwardArgs(args)...),
collator_error(U_ZERO_ERROR),
collator(Collator::createInstance(collator_error))
{
collator-setStrength(Collator::PRIMARY);
}
void node(const shared_ptrOsmium::OSM::Node n) {
do_tag_update(n-tags(), n-get_type());
Base::node(n);
}
void way(const shared_ptrOsmium::OSM::Way w) {
do_tag_update(w-tags(), w-get_type());
Base::way(w);
}
void relation(const shared_ptrOsmium::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
templatetypename Callback
class Multipolygon : public Base {
using std::vector;
/// a list of areas that need to be completed
vectorOsmium::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_maposm_object_id_t, vectorosm_object_id_t way2areaidx_t;
way2areaidx_t m_way2areaidx;
bool m_attempt_repair;
Callback
Re: [mkgmap-dev] Computing area
Oops - old version. Here's the correct one:
#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::SparseTableOsmium::OSM::Position storage_sparsetable_t;
typedef Osmium::Storage::MmapOsmium::OSM::Position storage_mmap_t;
typedef Osmium::Handler::CoordinatesForWaysstorage_sparsetable_t, storage_mmap_t cfw_handler_t;
templatetypename 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:
templatetypename... Args
ComputeAreaHandler(Args... args)
: Base(std::forwardArgs(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_ptrOsmium::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_ptrOsmium::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_ptrOsmium::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)))
+
