strategy/geographic/area_box.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2021, Oracle and/or its affiliates.
0004
0005 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0006
0007 // Licensed under the Boost Software License version 1.0.
0008 // http://www.boost.org/users/license.html
0009
0010 #ifndef BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_BOX_HPP
0011 #define BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_BOX_HPP
0012
0013
0014 #include <boost/geometry/core/radian_access.hpp>
0015 #include <boost/geometry/srs/spheroid.hpp>
0016 #include <boost/geometry/strategies/spherical/get_radius.hpp>
0017 #include <boost/geometry/strategy/area.hpp>
0018 #include <boost/geometry/util/normalize_spheroidal_box_coordinates.hpp>
0019
0020
0021 namespace boost { namespace geometry
0022 {
0023
0024 namespace strategy { namespace area
0025 {
0026
0027 // Based on the approach for spherical coordinate system:
0028 // https://math.stackexchange.com/questions/131735/surface-element-in-spherical-coordinates
0029 // http://www.cs.cmu.edu/afs/cs/academic/class/16823-s16/www/pdfs/appearance-modeling-3.pdf
0030 // https://www.astronomyclub.xyz/celestial-sphere-2/solid-angle-on-the-celestial-sphere.html
0031 // https://mathworld.wolfram.com/SolidAngle.html
0032 // https://en.wikipedia.org/wiki/Spherical_coordinate_system
0033 // and equations for spheroid:
0034 // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
0035 // https://en.wikipedia.org/wiki/Meridian_arc
0036 // Note that the equations use geodetic latitudes so we do not have to convert them.
0037 // assume(y_max > y_min);
0038 // assume(x_max > x_min);
0039 // M: a*(1-e^2) / (1-e^2*sin(y)^2)^(3/2);
0040 // N: a / sqrt(1-e^2*sin(y)^2);
0041 // O: N*cos(y)*M;
0042 // tellsimp(log(abs(e*sin(y_min)+1)), p_min);
0043 // tellsimp(log(abs(e*sin(y_min)-1)), m_min);
0044 // tellsimp(log(abs(e*sin(y_max)+1)), p_max);
0045 // tellsimp(log(abs(e*sin(y_max)-1)), m_max);
0046 // S: integrate(integrate(O, y, y_min, y_max), x, x_min, x_max);
0047 // combine(S);
0048 //
0049 // An alternative solution to the above formula was suggested by Charles Karney
0050 // https://github.com/boostorg/geometry/pull/832
0051 // The following are formulas for area of a box defined by the equator and some latitude,
0052 // not arbitrary box.
0053 // For e^2 > 0
0054 // dlambda*b^2*sin(phi)/2*(1/(1-e^2*sin(phi)^2) + atanh(e*sin(phi))/(e*sin(phi)))
0055 // For e^2 < 0
0056 // dlambda*b^2*sin(phi)/2*(1/(1-e^2*sin(phi)^2) + atan(ea*sin(phi))/(ea*sin(phi)))
0057 // where ea = sqrt(-e^2)
0058 template
0059 <
0060     typename Spheroid = srs::spheroid<double>,
0061     typename CalculationType = void
0062 >
0063 class geographic_box
0064 {
0065 public:
0066     template <typename Box>
0067     struct result_type
0068         : strategy::area::detail::result_type
0069             <
0070                 Box,
0071                 CalculationType
0072             >
0073     {};
0074
0075     geographic_box() = default;
0076
0077     explicit geographic_box(Spheroid const& spheroid)
0078         : m_spheroid(spheroid)
0079     {}
0080
0081     template <typename Box>
0082     inline auto apply(Box const& box) const
0083     {
0084         typedef typename result_type<Box>::type return_type;
0085
0086         return_type const c0 = 0;
0087
0088         return_type x_min = get_as_radian<min_corner, 0>(box); // lon
0089         return_type y_min = get_as_radian<min_corner, 1>(box); // lat
0090         return_type x_max = get_as_radian<max_corner, 0>(box);
0091         return_type y_max = get_as_radian<max_corner, 1>(box);
0092
0093         math::normalize_spheroidal_box_coordinates<radian>(x_min, y_min, x_max, y_max);
0094
0095         if (x_min == x_max || y_max == y_min)
0096         {
0097             return c0;
0098         }
0099
0100         return_type const e2 = formula::eccentricity_sqr<return_type>(m_spheroid);
0101
0102         return_type const x_diff = x_max - x_min;
0103         return_type const sin_y_min = sin(y_min);
0104         return_type const sin_y_max = sin(y_max);
0105
0106         if (math::equals(e2, c0))
0107         {
0108             // spherical formula
0109             return_type const a = get_radius<0>(m_spheroid);
0110             return x_diff * (sin_y_max - sin_y_min) * a * a;
0111         }
0112
0113         return_type const c1 = 1;
0114         return_type const c2 = 2;
0115         return_type const b = get_radius<2>(m_spheroid);
0116
0117         /*
0118         return_type const c4 = 4;
0119         return_type const e = math::sqrt(e2);
0120
0121         return_type const p_min = log(math::abs(e * sin_y_min + c1));
0122         return_type const p_max = log(math::abs(e * sin_y_max + c1));
0123         return_type const m_min = log(math::abs(e * sin_y_min - c1));
0124         return_type const m_max = log(math::abs(e * sin_y_max - c1));
0125         return_type const n_min = e * sin_y_min * sin_y_min;
0126         return_type const n_max = e * sin_y_max * sin_y_max;
0127         return_type const d_min = e * n_min - c1;
0128         return_type const d_max = e * n_max - c1;
0129
0130         // NOTE: For equal latitudes the original formula generated by maxima may give negative
0131         //   result. It's caused by the order of operations, so here they're rearranged for
0132         //   symmetry.
0133         return_type const comp0 = (p_min - m_min) / (c4 * e * d_min);
0134         return_type const comp1 = sin_y_min / (c2 * d_min);
0135         return_type const comp2 = n_min * (m_min - p_min) / (c4 * d_min);
0136         return_type const comp3 = (p_max - m_max) / (c4 * e * d_max);
0137         return_type const comp4 = sin_y_max / (c2 * d_max);
0138         return_type const comp5 = n_max * (m_max - p_max) / (c4 * d_max);
0139         return_type const comp02 = comp0 + comp1 + comp2;
0140         return_type const comp35 = comp3 + comp4 + comp5;
0141
0142         return b * b * x_diff * (comp02 - comp35);
0143         */
0144
0145         return_type const comp0_min = c1 / (c1 - e2 * sin_y_min * sin_y_min);
0146         return_type const comp0_max = c1 / (c1 - e2 * sin_y_max * sin_y_max);
0147
0148         // NOTE: For latitudes equal to 0 the original formula returns NAN
0149         return_type comp1_min = 0, comp1_max = 0;
0150         if (e2 > c0)
0151         {
0152             return_type const e = math::sqrt(e2);
0153             return_type const e_sin_y_min = e * sin_y_min;
0154             return_type const e_sin_y_max = e * sin_y_max;
0155
0156             comp1_min = e_sin_y_min == c0 ? c1 : atanh(e_sin_y_min) / e_sin_y_min;
0157             comp1_max = e_sin_y_max == c0 ? c1 : atanh(e_sin_y_max) / e_sin_y_max;
0158         }
0159         else
0160         {
0161             return_type const ea = math::sqrt(-e2);
0162             return_type const ea_sin_y_min = ea * sin_y_min;
0163             return_type const ea_sin_y_max = ea * sin_y_max;
0164
0165             comp1_min = ea_sin_y_min == c0 ? c1 : atan(ea_sin_y_min) / ea_sin_y_min;
0166             comp1_max = ea_sin_y_max == c0 ? c1 : atan(ea_sin_y_max) / ea_sin_y_max;
0167         }
0168
0169         return_type const comp01_min = sin_y_min * (comp0_min + comp1_min);
0170         return_type const comp01_max = sin_y_max * (comp0_max + comp1_max);
0171
0172         return b * b * x_diff * (comp01_max - comp01_min) / c2;
0173     }
0174
0175     Spheroid model() const
0176     {
0177         return m_spheroid;
0178     }
0179
0180 private:
0181     Spheroid m_spheroid;
0182 };
0183
0184
0185 }} // namespace strategy::area
0186
0187
0188 }} // namespace boost::geometry
0189
0190
0191 #endif // BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_BOX_HPP