192 lines
6.6 KiB
C++
192 lines
6.6 KiB
C++
// Boost.Geometry
|
|
|
|
// Copyright (c) 2021, Oracle and/or its affiliates.
|
|
|
|
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
|
|
|
// Licensed under the Boost Software License version 1.0.
|
|
// http://www.boost.org/users/license.html
|
|
|
|
#ifndef BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_BOX_HPP
|
|
#define BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_BOX_HPP
|
|
|
|
|
|
#include <boost/geometry/core/radian_access.hpp>
|
|
#include <boost/geometry/srs/spheroid.hpp>
|
|
#include <boost/geometry/strategies/spherical/get_radius.hpp>
|
|
#include <boost/geometry/strategy/area.hpp>
|
|
#include <boost/geometry/util/normalize_spheroidal_box_coordinates.hpp>
|
|
|
|
|
|
namespace boost { namespace geometry
|
|
{
|
|
|
|
namespace strategy { namespace area
|
|
{
|
|
|
|
// Based on the approach for spherical coordinate system:
|
|
// https://math.stackexchange.com/questions/131735/surface-element-in-spherical-coordinates
|
|
// http://www.cs.cmu.edu/afs/cs/academic/class/16823-s16/www/pdfs/appearance-modeling-3.pdf
|
|
// https://www.astronomyclub.xyz/celestial-sphere-2/solid-angle-on-the-celestial-sphere.html
|
|
// https://mathworld.wolfram.com/SolidAngle.html
|
|
// https://en.wikipedia.org/wiki/Spherical_coordinate_system
|
|
// and equations for spheroid:
|
|
// https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
|
|
// https://en.wikipedia.org/wiki/Meridian_arc
|
|
// Note that the equations use geodetic latitudes so we do not have to convert them.
|
|
// assume(y_max > y_min);
|
|
// assume(x_max > x_min);
|
|
// M: a*(1-e^2) / (1-e^2*sin(y)^2)^(3/2);
|
|
// N: a / sqrt(1-e^2*sin(y)^2);
|
|
// O: N*cos(y)*M;
|
|
// tellsimp(log(abs(e*sin(y_min)+1)), p_min);
|
|
// tellsimp(log(abs(e*sin(y_min)-1)), m_min);
|
|
// tellsimp(log(abs(e*sin(y_max)+1)), p_max);
|
|
// tellsimp(log(abs(e*sin(y_max)-1)), m_max);
|
|
// S: integrate(integrate(O, y, y_min, y_max), x, x_min, x_max);
|
|
// combine(S);
|
|
//
|
|
// An alternative solution to the above formula was suggested by Charles Karney
|
|
// https://github.com/boostorg/geometry/pull/832
|
|
// The following are formulas for area of a box defined by the equator and some latitude,
|
|
// not arbitrary box.
|
|
// For e^2 > 0
|
|
// dlambda*b^2*sin(phi)/2*(1/(1-e^2*sin(phi)^2) + atanh(e*sin(phi))/(e*sin(phi)))
|
|
// For e^2 < 0
|
|
// dlambda*b^2*sin(phi)/2*(1/(1-e^2*sin(phi)^2) + atan(ea*sin(phi))/(ea*sin(phi)))
|
|
// where ea = sqrt(-e^2)
|
|
template
|
|
<
|
|
typename Spheroid = srs::spheroid<double>,
|
|
typename CalculationType = void
|
|
>
|
|
class geographic_box
|
|
{
|
|
public:
|
|
template <typename Box>
|
|
struct result_type
|
|
: strategy::area::detail::result_type
|
|
<
|
|
Box,
|
|
CalculationType
|
|
>
|
|
{};
|
|
|
|
geographic_box() = default;
|
|
|
|
explicit geographic_box(Spheroid const& spheroid)
|
|
: m_spheroid(spheroid)
|
|
{}
|
|
|
|
template <typename Box>
|
|
inline auto apply(Box const& box) const
|
|
{
|
|
typedef typename result_type<Box>::type return_type;
|
|
|
|
return_type const c0 = 0;
|
|
|
|
return_type x_min = get_as_radian<min_corner, 0>(box); // lon
|
|
return_type y_min = get_as_radian<min_corner, 1>(box); // lat
|
|
return_type x_max = get_as_radian<max_corner, 0>(box);
|
|
return_type y_max = get_as_radian<max_corner, 1>(box);
|
|
|
|
math::normalize_spheroidal_box_coordinates<radian>(x_min, y_min, x_max, y_max);
|
|
|
|
if (x_min == x_max || y_max == y_min)
|
|
{
|
|
return c0;
|
|
}
|
|
|
|
return_type const e2 = formula::eccentricity_sqr<return_type>(m_spheroid);
|
|
|
|
return_type const x_diff = x_max - x_min;
|
|
return_type const sin_y_min = sin(y_min);
|
|
return_type const sin_y_max = sin(y_max);
|
|
|
|
if (math::equals(e2, c0))
|
|
{
|
|
// spherical formula
|
|
return_type const a = get_radius<0>(m_spheroid);
|
|
return x_diff * (sin_y_max - sin_y_min) * a * a;
|
|
}
|
|
|
|
return_type const c1 = 1;
|
|
return_type const c2 = 2;
|
|
return_type const b = get_radius<2>(m_spheroid);
|
|
|
|
/*
|
|
return_type const c4 = 4;
|
|
return_type const e = math::sqrt(e2);
|
|
|
|
return_type const p_min = log(math::abs(e * sin_y_min + c1));
|
|
return_type const p_max = log(math::abs(e * sin_y_max + c1));
|
|
return_type const m_min = log(math::abs(e * sin_y_min - c1));
|
|
return_type const m_max = log(math::abs(e * sin_y_max - c1));
|
|
return_type const n_min = e * sin_y_min * sin_y_min;
|
|
return_type const n_max = e * sin_y_max * sin_y_max;
|
|
return_type const d_min = e * n_min - c1;
|
|
return_type const d_max = e * n_max - c1;
|
|
|
|
// NOTE: For equal latitudes the original formula generated by maxima may give negative
|
|
// result. It's caused by the order of operations, so here they're rearranged for
|
|
// symmetry.
|
|
return_type const comp0 = (p_min - m_min) / (c4 * e * d_min);
|
|
return_type const comp1 = sin_y_min / (c2 * d_min);
|
|
return_type const comp2 = n_min * (m_min - p_min) / (c4 * d_min);
|
|
return_type const comp3 = (p_max - m_max) / (c4 * e * d_max);
|
|
return_type const comp4 = sin_y_max / (c2 * d_max);
|
|
return_type const comp5 = n_max * (m_max - p_max) / (c4 * d_max);
|
|
return_type const comp02 = comp0 + comp1 + comp2;
|
|
return_type const comp35 = comp3 + comp4 + comp5;
|
|
|
|
return b * b * x_diff * (comp02 - comp35);
|
|
*/
|
|
|
|
return_type const comp0_min = c1 / (c1 - e2 * sin_y_min * sin_y_min);
|
|
return_type const comp0_max = c1 / (c1 - e2 * sin_y_max * sin_y_max);
|
|
|
|
// NOTE: For latitudes equal to 0 the original formula returns NAN
|
|
return_type comp1_min = 0, comp1_max = 0;
|
|
if (e2 > c0)
|
|
{
|
|
return_type const e = math::sqrt(e2);
|
|
return_type const e_sin_y_min = e * sin_y_min;
|
|
return_type const e_sin_y_max = e * sin_y_max;
|
|
|
|
comp1_min = e_sin_y_min == c0 ? c1 : atanh(e_sin_y_min) / e_sin_y_min;
|
|
comp1_max = e_sin_y_max == c0 ? c1 : atanh(e_sin_y_max) / e_sin_y_max;
|
|
}
|
|
else
|
|
{
|
|
return_type const ea = math::sqrt(-e2);
|
|
return_type const ea_sin_y_min = ea * sin_y_min;
|
|
return_type const ea_sin_y_max = ea * sin_y_max;
|
|
|
|
comp1_min = ea_sin_y_min == c0 ? c1 : atan(ea_sin_y_min) / ea_sin_y_min;
|
|
comp1_max = ea_sin_y_max == c0 ? c1 : atan(ea_sin_y_max) / ea_sin_y_max;
|
|
}
|
|
|
|
return_type const comp01_min = sin_y_min * (comp0_min + comp1_min);
|
|
return_type const comp01_max = sin_y_max * (comp0_max + comp1_max);
|
|
|
|
return b * b * x_diff * (comp01_max - comp01_min) / c2;
|
|
}
|
|
|
|
Spheroid model() const
|
|
{
|
|
return m_spheroid;
|
|
}
|
|
|
|
private:
|
|
Spheroid m_spheroid;
|
|
};
|
|
|
|
|
|
}} // namespace strategy::area
|
|
|
|
|
|
}} // namespace boost::geometry
|
|
|
|
|
|
#endif // BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_BOX_HPP
|