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::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;
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);
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;
}
Sorry!
_______________________________________________
mkgmap-dev mailing list
mkgmap-dev@lists.mkgmap.org.uk
http://lists.mkgmap.org.uk/mailman/listinfo/mkgmap-dev
