1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
|
// 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
|
