From 3fda39ded1a3d89a7eeeddb1e0ff786c6684a08b Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 8 Nov 2016 13:02:37 +0200 Subject: [PATCH] [vertex latitude] Simplify formulas --- .../geometry/formulas/vertex_latitude.hpp | 27 ++++ .../envelope_expand/envelope_on_spheroid.cpp | 126 +++++++++--------- 2 files changed, 90 insertions(+), 63 deletions(-) diff --git a/include/boost/geometry/formulas/vertex_latitude.hpp b/include/boost/geometry/formulas/vertex_latitude.hpp index f415d7c42..d55e15c4a 100644 --- a/include/boost/geometry/formulas/vertex_latitude.hpp +++ b/include/boost/geometry/formulas/vertex_latitude.hpp @@ -53,6 +53,11 @@ public: typename T1, typename T2 > +/* + * formula based on paper + * [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4), + * 637–644, 1996 + static inline CT apply(T1 const& lat1, T2 const& alp1) { @@ -71,6 +76,28 @@ public: / (e2_sin2 - e2 * cos2_sin2))); return vertex_lat; } +*/ + + // simpler formula based on Clairaut relation for spheroids + + static inline CT apply(T1 const& lat1, + T2 const& alp1) + { + + geometry::srs::spheroid spheroid; + CT const f = detail::flattening(spheroid); + + CT const one_minus_f = (CT(1) - f); + + //get the reduced latitude + CT const bet1 = atan( one_minus_f * tan(lat1) ); + + //apply Clairaut relation + CT const betv = vertex_latitude_on_sphere::apply(bet1, alp1); + + //return the spheroid latitude + return atan( tan(betv) / one_minus_f ); + } /* template diff --git a/test/algorithms/envelope_expand/envelope_on_spheroid.cpp b/test/algorithms/envelope_expand/envelope_on_spheroid.cpp index c9bca08fc..a48562d6f 100644 --- a/test/algorithms/envelope_expand/envelope_on_spheroid.cpp +++ b/test/algorithms/envelope_expand/envelope_on_spheroid.cpp @@ -15,7 +15,7 @@ #include #ifndef BOOST_TEST_MODULE -#define BOOST_TEST_MODULE test_envelope_on_spheroid +#define BOOST_TEST_MODULE test_envelope_on_sphere_or_spheroid #endif #include @@ -288,7 +288,7 @@ template template class Inverse = bg::formula::thomas_inverse > -struct test_envelope_on_spheroid +struct test_envelope_on_sphere_or_spheroid { static inline void apply(std::string const& case_id, Geometry const& geometry, @@ -318,7 +318,7 @@ struct test_envelope_on_spheroid >::apply(reversed_case_id, reversed_geometry, lon_min2, lat_min2, height_min2, lon_max2, lat_max2, height_max2, - tolerance); + tolerance); } #ifdef BOOST_GEOMETRY_TEST_DEBUG @@ -370,7 +370,7 @@ struct test_envelope_on_spheroid // special tester for rings template -struct test_envelope_on_spheroid +struct test_envelope_on_sphere_or_spheroid { static inline void apply(std::string const& case_id, Geometry const& geometry, @@ -394,7 +394,7 @@ struct test_envelope_on_spheroid typename bg::point_type::type, false > ccw_ring; bg::convert(geometry, ccw_ring); - + envelope_on_spheroid_basic_tester < MBR @@ -429,7 +429,7 @@ void test_empty_geometry(std::string const& case_id, std::string const& wkt) typedef bg::model::point point_type; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; typedef typename bg::coordinate_type::type ct; ct high_val = boost::numeric::bounds::highest(); @@ -456,7 +456,7 @@ void test_envelope_point() typedef bg::model::point point_type; typedef point_type G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; tester::apply("p01", from_wkt("POINT(10 10)"), @@ -543,7 +543,7 @@ void test_envelope_point_with_height() typedef bg::model::point point_type; typedef point_type G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; tester::apply("ph01", from_wkt("POINT(10 10 1256)"), @@ -566,7 +566,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_sphere ) typedef bg::model::point P; typedef bg::model::segment

G; typedef bg::model::box

B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; double const eps = std::numeric_limits::epsilon(); @@ -576,7 +576,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_sphere ) tester::apply("s02", from_wkt("SEGMENT(10 10,40 10)"), - 10, 10, 40, 10.34527004614999); + 10, 10, 40, 10.345270046149988); tester::apply("s02a", from_wkt("SEGMENT(40 10,10 10)"), @@ -748,7 +748,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid ) typedef bg::model::point P; typedef bg::model::segment

G; typedef bg::model::box

B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; double const eps = std::numeric_limits::epsilon(); @@ -758,23 +758,23 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid ) tester::apply("s02", from_wkt("SEGMENT(10 10,40 10)"), - 10, 10, 40, 10.347587605817935); + 10, 10, 40, 10.347587605817942); tester::apply("s02a", from_wkt("SEGMENT(40 10,10 10)"), - 10, 10, 40, 10.347587605817935); + 10, 10, 40, 10.347587605817942); tester::apply("s03", from_wkt("SEGMENT(160 10,-170 10)"), - 160, 10, 190, 10.347587605817935); + 160, 10, 190, 10.347587605817942); tester::apply("s03a", from_wkt("SEGMENT(-170 10,160 10)"), - 160, 10, 190, 10.347587605817935); + 160, 10, 190, 10.347587605817942); tester::apply("s03b", from_wkt("SEGMENT(-170 -10,160 -10)"), - 160, -10.347587605817935, 190, -10); + 160, -10.347587605817942, 190, -10); tester::apply("s04", from_wkt("SEGMENT(-40 45,140 60)"), @@ -826,12 +826,12 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid ) // very long segment tester::apply("s10", from_wkt("SEGMENT(0 -45,181 30)"), - -179, -88.111912243387664, 0, 30, + -179, -88.111912243387607, 0, 30, 2.0 * eps); tester::apply("s11", from_wkt("SEGMENT(260 30,20 45)"), - -100, 30, 20, 57.990810958016951); + -100, 30, 20, 57.990810958016965); tester::apply("s11a", from_wkt("SEGMENT(260 45,20 30)"), @@ -918,7 +918,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid ) 1-heps, 1, 1, 1); tester::apply("s104", G(P(2, 1), P(1, 1-heps)), - 1, 1-heps, 2, 1.000038327157039); + 1, 1-heps, 2, 1.0000383271569036); tester::apply("s105", G(P(1, 2), P(1-heps, 1)), 1-heps, 1, 1, 2); @@ -931,7 +931,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_thomas ) typedef bg::model::point P; typedef bg::model::segment

G; typedef bg::model::box

B; - typedef test_envelope_on_spheroid + typedef test_envelope_on_sphere_or_spheroid < G, B, typename bg::tag::type, @@ -947,23 +947,23 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_thomas ) tester::apply("s02", from_wkt("SEGMENT(10 10,40 10)"), - 10, 10, 40, 10.347587605817935); + 10, 10, 40, 10.347587605817942); tester::apply("s02a", from_wkt("SEGMENT(40 10,10 10)"), - 10, 10, 40, 10.347587605817935); + 10, 10, 40, 10.347587605817942); tester::apply("s03", from_wkt("SEGMENT(160 10,-170 10)"), - 160, 10, 190, 10.347587605817935); + 160, 10, 190, 10.347587605817942); tester::apply("s03a", from_wkt("SEGMENT(-170 10,160 10)"), - 160, 10, 190, 10.347587605817935); + 160, 10, 190, 10.347587605817942); tester::apply("s03b", from_wkt("SEGMENT(-170 -10,160 -10)"), - 160, -10.347587605817935, 190, -10); + 160, -10.347587605817942, 190, -10); tester::apply("s04", from_wkt("SEGMENT(-40 45,140 60)"), @@ -1015,12 +1015,12 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_thomas ) // very long segment tester::apply("s10", from_wkt("SEGMENT(0 -45,181 30)"), - -179, -88.111912243387664, 0, 30, + -179, -88.111912243387607, 0, 30, 2.0 * eps); tester::apply("s11", from_wkt("SEGMENT(260 30,20 45)"), - -100, 30, 20, 57.990810958016951); + -100, 30, 20, 57.990810958016965); tester::apply("s11a", from_wkt("SEGMENT(260 45,20 30)"), @@ -1089,7 +1089,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_andoyer ) typedef bg::model::point P; typedef bg::model::segment

G; typedef bg::model::box

B; - typedef test_envelope_on_spheroid + typedef test_envelope_on_sphere_or_spheroid < G, B, typename bg::tag::type, @@ -1105,23 +1105,23 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_andoyer ) tester::apply("s02", from_wkt("SEGMENT(10 10,40 10)"), - 10, 10, 40, 10.347587099602052); + 10, 10, 40, 10.34758709960203); tester::apply("s02a", from_wkt("SEGMENT(40 10,10 10)"), - 10, 10, 40, 10.347587099602052); + 10, 10, 40, 10.34758709960203); tester::apply("s03", from_wkt("SEGMENT(160 10,-170 10)"), - 160, 10, 190, 10.347587099602052); + 160, 10, 190, 10.34758709960203); tester::apply("s03a", from_wkt("SEGMENT(-170 10,160 10)"), - 160, 10, 190, 10.347587099602052); + 160, 10, 190, 10.34758709960203); tester::apply("s03b", from_wkt("SEGMENT(-170 -10,160 -10)"), - 160, -10.347587099602052, 190, -10); + 160, -10.34758709960203, 190, -10); tester::apply("s04", from_wkt("SEGMENT(-40 45,140 60)"), @@ -1168,7 +1168,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_andoyer ) tester::apply("s09a", from_wkt("SEGMENT(2 -45,181 30)"), - 2, -87.690317839849669, 181, 30); + 2, -87.690317839849726, 181, 30); // very long segment tester::apply("s10", @@ -1178,7 +1178,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_andoyer ) tester::apply("s11", from_wkt("SEGMENT(260 30,20 45)"), - -100, 30, 20, 57.990742552280139); + -100, 30, 20, 57.990742552280153); tester::apply("s11a", from_wkt("SEGMENT(260 45,20 30)"), @@ -1247,7 +1247,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_vincenty ) typedef bg::model::point P; typedef bg::model::segment

G; typedef bg::model::box

B; - typedef test_envelope_on_spheroid + typedef test_envelope_on_sphere_or_spheroid < G, B, typename bg::tag::type, @@ -1263,23 +1263,23 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_vincenty ) tester::apply("s02", from_wkt("SEGMENT(10 10,40 10)"), - 10, 10, 40, 10.347587628821922); + 10, 10, 40, 10.347587628821941); tester::apply("s02a", from_wkt("SEGMENT(40 10,10 10)"), - 10, 10, 40, 10.347587628821922); + 10, 10, 40, 10.347587628821941); tester::apply("s03", from_wkt("SEGMENT(160 10,-170 10)"), - 160, 10, 190, 10.347587628821922); + 160, 10, 190, 10.347587628821941); tester::apply("s03a", from_wkt("SEGMENT(-170 10,160 10)"), - 160, 10, 190, 10.347587628821922); + 160, 10, 190, 10.347587628821941); tester::apply("s03b", from_wkt("SEGMENT(-170 -10,160 -10)"), - 160, -10.347587628821922, 190, -10); + 160, -10.347587628821941, 190, -10); tester::apply("s04", from_wkt("SEGMENT(-40 45,140 60)"), @@ -1331,7 +1331,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_vincenty ) // very long segment tester::apply("s10", from_wkt("SEGMENT(0 -45,181 30)"), - -179, -88.11193623291669, 0, 30, + -179, -88.111936232916634, 0, 30, 2.0 * eps); tester::apply("s11", @@ -1358,7 +1358,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_strategy_vincenty ) tester::apply("s15", from_wkt("SEGMENT(50 45,185 45)"), - 50, 45, 185, 69.098479136978511); + 50, 45, 185, 69.098479136978497); // segment that lies on the equator tester::apply("s16", @@ -1404,7 +1404,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_sphere_with_height ) typedef bg::model::point point_type; typedef bg::model::segment G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; tester::apply("sh01", from_wkt("SEGMENT(10 10 567,40 40 1356)"), @@ -1421,7 +1421,7 @@ BOOST_AUTO_TEST_CASE( envelope_segment_spheroid_with_height ) typedef bg::model::point point_type; typedef bg::model::segment G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; tester::apply("sh01", from_wkt("SEGMENT(10 10 567,40 40 1356)"), @@ -1438,7 +1438,7 @@ void test_envelope_multipoint() typedef bg::model::point P; typedef bg::model::multi_point

G; typedef bg::model::box

B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; // empty multipoint test_empty_geometry("mp00", "MULTIPOINT()"); @@ -1590,7 +1590,7 @@ void test_envelope_multipoint_with_height() typedef bg::model::point point_type; typedef bg::model::multi_point G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; // empty multipoint test_empty_geometry("mph00", "MULTIPOINT()"); @@ -1621,7 +1621,7 @@ void test_envelope_box() typedef bg::model::point P; typedef bg::model::box

G; typedef bg::model::box

B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; tester::apply("b01", from_wkt("BOX(10 10,20 20)"), @@ -1821,7 +1821,7 @@ void test_envelope_box_with_height() typedef bg::model::point point_type; typedef bg::model::box G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; tester::apply("bh01", from_wkt("BOX(10 10 567,20 20 2834)"), @@ -1849,7 +1849,7 @@ BOOST_AUTO_TEST_CASE( envelope_sphere_linestring ) typedef bg::model::point P; typedef bg::model::linestring

G; typedef bg::model::box

B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; // empty linestring test_empty_geometry("l00", "LINESTRING()"); @@ -1982,7 +1982,7 @@ BOOST_AUTO_TEST_CASE( envelope_spheroid_linestring ) typedef bg::model::point P; typedef bg::model::linestring

G; typedef bg::model::box

B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; // empty linestring test_empty_geometry("l00", "LINESTRING()"); @@ -2047,7 +2047,7 @@ BOOST_AUTO_TEST_CASE( envelope_spheroid_linestring ) // linestring that crosses the antimeridian tester::apply("l06a", from_wkt("LINESTRING(-130 85,-170 84,170 40,160 80)"), - 160, 40, 230, 85.026305563151482); + 160, 40, 230, 85.02630556315151); // linestring that goes through the north pole tester::apply("l07", @@ -2116,7 +2116,7 @@ BOOST_AUTO_TEST_CASE( envelope_linestring_sphere_with_height ) typedef bg::model::point point_type; typedef bg::model::linestring G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; // empty linestring test_empty_geometry("lh00", "LINESTRING()"); @@ -2132,7 +2132,7 @@ BOOST_AUTO_TEST_CASE( envelope_linestring_spheroid_with_height ) typedef bg::model::point point_type; typedef bg::model::linestring G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; // empty linestring test_empty_geometry("lh00", "LINESTRING()"); @@ -2148,7 +2148,7 @@ BOOST_AUTO_TEST_CASE( envelope_sphere_multilinestring ) typedef bg::model::point point_type; typedef bg::model::multi_linestring > G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; // empty multilinestring test_empty_geometry("ml00", "MULTILINESTRING()"); @@ -2212,7 +2212,7 @@ BOOST_AUTO_TEST_CASE( envelope_spheroid_multilinestring ) typedef bg::model::point point_type; typedef bg::model::multi_linestring > G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; // empty multilinestring test_empty_geometry("ml00", "MULTILINESTRING()"); @@ -2277,7 +2277,7 @@ BOOST_AUTO_TEST_CASE( envelope_multilinestring_sphere_with_height ) typedef bg::model::point point_type; typedef bg::model::multi_linestring > G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; tester::apply("mlh01", from_wkt("MULTILINESTRING((10 15 1000))"), @@ -2300,7 +2300,7 @@ BOOST_AUTO_TEST_CASE( envelope_multilinestring_spheroid_with_height ) typedef bg::model::point point_type; typedef bg::model::multi_linestring > G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; tester::apply("mlh01", from_wkt("MULTILINESTRING((10 15 1000))"), @@ -2318,8 +2318,8 @@ BOOST_AUTO_TEST_CASE( envelope_multilinestring_spheroid_with_height ) } -// unit test for rings that cover more than half of the earth is de-activated -// for now since these rings are not supported +// unit test for rings de-activated for now (current implementation +// for area on the spherical equatorial coordinate system is not complete) // TODO: re-activate once implementation is done BOOST_AUTO_TEST_CASE( envelope_cw_ring ) { @@ -2327,14 +2327,14 @@ BOOST_AUTO_TEST_CASE( envelope_cw_ring ) typedef bg::model::point point_type; typedef bg::model::polygon G; typedef bg::model::box B; - typedef test_envelope_on_spheroid tester; + typedef test_envelope_on_sphere_or_spheroid tester; + + //double const eps = std::numeric_limits::epsilon(); tester::apply("r01cw", from_wkt("POLYGON((0 10,0 45,50 10,0 10))"), 0, 10, 50, 45); #if 0 - double const eps = std::numeric_limits::epsilon(); - // ring that contains both the north and south poles in its interior tester::apply("r01cw-r", from_wkt("POLYGON((0 10,50 10,0 45,0 10))"),