mirror of
https://github.com/boostorg/geometry.git
synced 2026-01-27 06:52:11 +00:00
184 lines
6.3 KiB
C++
184 lines
6.3 KiB
C++
// Boost.Geometry (aka GGL, Generic Geometry Library)
|
|
// Unit Test
|
|
|
|
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
|
|
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
|
|
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
|
|
|
|
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
|
|
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
|
|
|
|
// Use, modification and distribution is subject to the Boost Software License,
|
|
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
|
// http://www.boost.org/LICENSE_1_0.txt)
|
|
|
|
#include <geometry_test_common.hpp>
|
|
|
|
#include <boost/geometry/algorithms/area.hpp>
|
|
#include <boost/geometry/algorithms/transform.hpp>
|
|
#include <boost/geometry/geometries/geometries.hpp>
|
|
#include <boost/geometry/geometries/point_xy.hpp>
|
|
#include <boost/geometry/io/wkt/read.hpp>
|
|
#include <boost/geometry/strategies/strategies.hpp>
|
|
|
|
#include <boost/geometry/extensions/gis/latlong/point_ll.hpp>
|
|
#include <boost/geometry/extensions/gis/geographic/strategies/area_huiller_earth.hpp>
|
|
#include <boost/geometry/extensions/gis/projections/epsg.hpp>
|
|
#include <boost/geometry/extensions/gis/projections/parameters.hpp>
|
|
#include <boost/geometry/extensions/gis/projections/proj/sterea.hpp>
|
|
|
|
//#include <test_common/test_point.hpp>
|
|
|
|
template <typename PRJ, typename XY, typename LL>
|
|
void add_to_ring(PRJ const& prj, LL const& ll,
|
|
bg::model::ring<LL>& ring_ll,
|
|
bg::model::ring<XY>& ring_xy)
|
|
{
|
|
ring_ll.push_back(ll);
|
|
|
|
XY xy;
|
|
prj.forward(ll, xy);
|
|
ring_xy.push_back(xy);
|
|
}
|
|
|
|
template <typename XY, typename LL>
|
|
void test_area_polygon_ll(bool concave, bool hole, double perc)
|
|
{
|
|
BOOST_ASSERT(! (concave && hole) );
|
|
|
|
typedef typename bg::coordinate_type<LL>::type T;
|
|
|
|
// Amsterdam, Rotterdam, The Hague, Utrecht,
|
|
// these cities together are the city group "Randstad"
|
|
|
|
LL a, r, h, u;
|
|
// Amsterdam 52 22'23"N 4 53'32"E
|
|
a.lat(bg::dms<bg::north, T>(52, 22, 23));
|
|
a.lon(bg::dms<bg::east, T>(4, 53, 32));
|
|
|
|
// Rotterdam 51 55'51"N 4 28'45"E
|
|
r.lat(bg::dms<bg::north, T>(51, 55, 51));
|
|
r.lon(bg::dms<bg::east, T>(4, 28, 45));
|
|
|
|
// The hague: 52 4' 48" N, 4 18' 0" E
|
|
h.lat(bg::dms<bg::north, T>(52, 4, 48));
|
|
h.lon(bg::dms<bg::east, T>(4, 18, 0));
|
|
|
|
// Utrecht
|
|
u.lat(bg::dms<bg::north, T>(52, 5, 36));
|
|
u.lon(bg::dms<bg::east, T>(5, 7, 10));
|
|
|
|
|
|
// For checking calculated area, use the Dutch projection (RD), this is EPSG code 28992
|
|
bg::projections::sterea_ellipsoid<LL, XY> dutch_prj(bg::projections::init(28992));
|
|
|
|
// Add them in clockwise direction
|
|
bg::model::polygon<LL> randstad;
|
|
bg::model::polygon<XY> randstad_xy;
|
|
add_to_ring(dutch_prj, a, randstad.outer(), randstad_xy.outer());
|
|
add_to_ring(dutch_prj, u, randstad.outer(), randstad_xy.outer());
|
|
|
|
// Concave case
|
|
if (concave)
|
|
{
|
|
// Add the city "Alphen" to create a concave case
|
|
// Alphen 52 7' 48" N, 4 39' 0" E
|
|
LL alphen(
|
|
bg::latitude<T>(bg::dms<bg::north, T>(52, 7, 48)),
|
|
bg::longitude<T>(bg::dms<bg::east, T>(4, 39)));
|
|
add_to_ring(dutch_prj, alphen, randstad.outer(), randstad_xy.outer());
|
|
}
|
|
|
|
add_to_ring(dutch_prj, r, randstad.outer(), randstad_xy.outer());
|
|
add_to_ring(dutch_prj, h, randstad.outer(), randstad_xy.outer());
|
|
add_to_ring(dutch_prj, a, randstad.outer(), randstad_xy.outer());
|
|
|
|
// Hole case
|
|
if (hole)
|
|
{
|
|
// Gouda 52 1' 12" N, 4 42' 0" E
|
|
LL gouda(
|
|
bg::latitude<T>(bg::dms<bg::north, T>(52, 1, 12)),
|
|
bg::longitude<T>(bg::dms<bg::east, T>(4, 42)));
|
|
// Woerden 52 5' 9" N, 4 53' 0" E
|
|
LL woerden(
|
|
bg::latitude<T>(bg::dms<bg::north, T>(52, 5, 9)),
|
|
bg::longitude<T>(bg::dms<bg::east, T>(4, 53, 0)));
|
|
// Uithoorn 52 13' 48" N, 4 49' 48" E
|
|
LL uithoorn(bg::latitude<T>
|
|
(bg::dms<bg::north, T>(52, 13, 48)),
|
|
bg::longitude<T>(bg::dms<bg::east, T>(4, 49, 48)));
|
|
// Alphen 52 7' 48" N, 4 39' 0" E
|
|
LL alphen(bg::latitude<T>(
|
|
bg::dms<bg::north, T>(52, 7, 48)),
|
|
bg::longitude<T>(bg::dms<bg::east, T>(4, 39)));
|
|
|
|
randstad.inners().resize(1);
|
|
randstad_xy.inners().resize(1);
|
|
|
|
typename bg::model::polygon<LL>::ring_type& ring = randstad.inners()[0];
|
|
typename bg::model::polygon<XY>::ring_type& ring_xy = randstad_xy.inners()[0];
|
|
|
|
// Add them in counter-clockwise direction (see map of the Netherlands)
|
|
add_to_ring(dutch_prj, gouda, ring, ring_xy);
|
|
add_to_ring(dutch_prj, woerden, ring, ring_xy);
|
|
add_to_ring(dutch_prj, uithoorn, ring, ring_xy);
|
|
add_to_ring(dutch_prj, alphen, ring, ring_xy);
|
|
add_to_ring(dutch_prj, gouda, ring, ring_xy);
|
|
}
|
|
|
|
|
|
// Check the area in square KM
|
|
static const double KM2 = 1.0e6;
|
|
double d_ll = bg::area(randstad) / KM2;
|
|
double d_xy = bg::area(randstad_xy) / KM2;
|
|
|
|
BOOST_CHECK_CLOSE(d_ll, d_xy, 1.0);
|
|
if (hole)
|
|
{
|
|
BOOST_CHECK_CLOSE(d_ll, 1148.210, perc);
|
|
BOOST_CHECK_CLOSE(d_xy, 1151.573, perc);
|
|
}
|
|
else
|
|
{
|
|
BOOST_CHECK_CLOSE(d_ll, concave ? 977.786 : 1356.168, perc);
|
|
BOOST_CHECK_CLOSE(d_xy, concave ? 980.658 : 1360.140, perc);
|
|
|
|
// No hole: area of outer should be equal to area of ring
|
|
double r_ll = bg::area(randstad.outer()) / KM2;
|
|
double r_xy = bg::area(randstad_xy.outer()) / KM2;
|
|
|
|
BOOST_CHECK_CLOSE(d_ll, r_ll, perc);
|
|
BOOST_CHECK_CLOSE(d_xy, r_xy, perc);
|
|
}
|
|
|
|
// Calculate are using specified strategy, here with radius in KM
|
|
// We then don't have to divide by KM*KM to get the same result
|
|
bg::strategy::area::huiller<LL, long double> strategy(6372.8);
|
|
d_ll = bg::area(randstad, strategy);
|
|
BOOST_CHECK_CLOSE(d_ll, d_xy, 1.0);
|
|
}
|
|
|
|
|
|
|
|
template <typename T>
|
|
void test_latlong(double perc)
|
|
{
|
|
test_area_polygon_ll<bg::model::d2::point_xy<T>, bg::model::ll::point<bg::degree, T> >(false, false, perc);
|
|
|
|
// with concavities
|
|
test_area_polygon_ll<bg::model::d2::point_xy<T>, bg::model::ll::point<bg::degree, T> >(true, false, perc);
|
|
|
|
// with holes
|
|
test_area_polygon_ll<bg::model::d2::point_xy<T>, bg::model::ll::point<bg::degree, T> >(false, true, perc);
|
|
|
|
}
|
|
|
|
int test_main(int, char* [])
|
|
{
|
|
test_latlong<double>(0.01);
|
|
//test_latlong<float>(0.3); // LL area calculations using projections differ
|
|
|
|
return 0;
|
|
}
|