GRASS GIS 7 Programmer's Manual
7.8.2(2019)-exported
geodist.c
Go to the documentation of this file.
1
/*!
2
* \file lib/gis/geodist.c
3
*
4
* \brief GIS Library - Geodesic distance routines.
5
*
6
* Distance from point to point along a geodesic code from Paul
7
* D. Thomas, 1970 "Spheroidal Geodesics, Reference Systems, and Local
8
* Geometry" U.S. Naval Oceanographic Office, p. 162 Engineering
9
* Library 526.3 T36s
10
* http://stinet.dtic.mil/oai/oai?&verb=getRecord&metadataPrefix=html&identifier=AD0703541
11
*
12
* <b>WARNING:</b> this code is preliminary and may be changed,
13
* including calling sequences to any of the functions defined here.
14
*
15
* (C) 2001-2009 by the GRASS Development Team
16
*
17
* This program is free software under the GNU General Public License
18
* (>=v2). Read the file COPYING that comes with GRASS for details.
19
*
20
* \author Original author CERL
21
*/
22
23
#include <math.h>
24
#include <grass/gis.h>
25
#include "
pi.h
"
26
27
static
struct
state
{
28
double
boa;
29
double
f;
30
double
ff64;
31
double
al;
32
double
t1, t2, t3, t4, t1r, t2r;
33
}
state
;
34
35
static
struct
state
*
st
= &
state
;
36
37
/*!
38
* \brief Begin geodesic distance.
39
*
40
* Initializes the distance calculations for the ellipsoid with
41
* semi-major axis <i>a</i> (in meters) and ellipsoid eccentricity squared
42
* <i>e2</i>. It is used only for the latitude-longitude projection.
43
*
44
* <b>Note:</b> Must be called once to establish the ellipsoid.
45
*
46
* \param a semi-major axis in meters
47
* \param e2 ellipsoid eccentricity
48
*/
49
50
void
G_begin_geodesic_distance
(
double
a,
double
e2)
51
{
52
st
->al = a;
53
st
->boa = sqrt(1 - e2);
54
st
->f = 1 -
st
->boa;
55
st
->ff64 =
st
->f *
st
->f / 64;
56
}
57
58
/*!
59
* \brief Sets geodesic distance lat1.
60
*
61
* Set the first latitude.
62
*
63
* <b>Note:</b> Must be called first.
64
*
65
* \param lat1 first latitude
66
* \return
67
*/
68
69
void
G_set_geodesic_distance_lat1
(
double
lat1)
70
{
71
st
->t1r = atan(
st
->boa * tan(
Radians
(lat1)));
72
}
73
74
/*!
75
* \brief Sets geodesic distance lat2.
76
*
77
* Set the second latitude.
78
*
79
* <b>Note:</b> Must be called second.
80
*
81
* \param lat2 second latitidue
82
*/
83
void
G_set_geodesic_distance_lat2
(
double
lat2)
84
{
85
double
stm, ctm, sdtm, cdtm;
86
double
tm, dtm;
87
88
st
->t2r = atan(
st
->boa * tan(
Radians
(lat2)));
89
90
tm = (
st
->t1r +
st
->t2r) / 2;
91
dtm = (
st
->t2r -
st
->t1r) / 2;
92
93
stm = sin(tm);
94
ctm = cos(tm);
95
sdtm = sin(dtm);
96
cdtm = cos(dtm);
97
98
st
->t1 = stm * cdtm;
99
st
->t1 =
st
->t1 *
st
->t1 * 2;
100
101
st
->t2 = sdtm * ctm;
102
st
->t2 =
st
->t2 *
st
->t2 * 2;
103
104
st
->t3 = sdtm * sdtm;
105
st
->t4 = cdtm * cdtm - stm * stm;
106
}
107
108
/*!
109
* \brief Calculates geodesic distance.
110
*
111
* Calculates the geodesic distance from <i>lon1,lat1</i> to
112
* <i>lon2,lat2</i> in meters where <i>lat1</i> was the latitude
113
* passed to G_set_geodesic_distance_latl() and <i>lat2</i> was the
114
* latitude passed to G_set_geodesic_distance_lat2().
115
*
116
* \param lon1 first longitude
117
* \param lon2 second longitude
118
*
119
* \return double distance in meters
120
*/
121
double
G_geodesic_distance_lon_to_lon
(
double
lon1,
double
lon2)
122
{
123
double
a, cd, d, e,
/*dl, */
124
q, sd, sdlmr,
t
, u, v,
x
, y;
125
126
127
sdlmr = sin(
Radians
(lon2 - lon1) / 2);
128
129
/* special case - shapiro */
130
if
(sdlmr == 0.0 &&
st
->t1r ==
st
->t2r)
131
return
0.0;
132
133
q =
st
->t3 + sdlmr * sdlmr *
st
->t4;
134
135
/* special case - shapiro */
136
if
(q == 1.0)
137
return
M_PI *
st
->al;
138
139
/* Mod: shapiro
140
* cd=1-2q is ill-conditioned if q is small O(10**-23)
141
* (for high lats? with lon1-lon2 < .25 degrees?)
142
* the computation of cd = 1-2*q will give cd==1.0.
143
* However, note that t=dl/sd is dl/sin(dl) which approaches 1 as dl->0.
144
* So the first step is to compute a good value for sd without using sin()
145
* and then check cd && q to see if we got cd==1.0 when we shouldn't.
146
* Note that dl isn't used except to get t,
147
* but both cd and sd are used later
148
*/
149
150
/* original code
151
cd=1-2*q;
152
dl=acos(cd);
153
sd=sin(dl);
154
t=dl/sd;
155
*/
156
157
cd = 1 - 2 * q;
/* ill-conditioned subtraction for small q */
158
/* mod starts here */
159
sd = 2 * sqrt(q - q * q);
/* sd^2 = 1 - cd^2 */
160
if
(q != 0.0 && cd == 1.0)
/* test for small q */
161
t
= 1.0;
162
else
if
(sd == 0.0)
163
t
= 1.0;
164
else
165
t
= acos(cd) / sd;
/* don't know how to fix acos(1-2*q) yet */
166
/* mod ends here */
167
168
u =
st
->t1 / (1 - q);
169
v =
st
->t2 / q;
170
d = 4 *
t
*
t
;
171
x
= u + v;
172
e = -2 * cd;
173
y = u - v;
174
a = -d * e;
175
176
return
st
->al * sd * (
t
177
-
st
->f / 4 * (
t
*
x
- y)
178
+
st
->ff64 * (
x
* (a + (
t
- (a + e) / 2) *
x
)
179
+ y * (-2 * d + e * y) + d *
x
* y));
180
}
181
182
/*!
183
* \brief Calculates geodesic distance.
184
*
185
* Calculates the geodesic distance from <i>lon1,lat1</i> to
186
* <i>lon2,lat2</i> in meters.
187
*
188
* <b>Note:</b> The calculation of the geodesic distance is fairly
189
* costly.
190
*
191
* \param lon1,lat1 longitude,latitude of first point
192
* \param lon2,lat2 longitude,latitude of second point
193
*
194
* \return distance in meters
195
*/
196
double
G_geodesic_distance
(
double
lon1,
double
lat1,
double
lon2,
double
lat2)
197
{
198
G_set_geodesic_distance_lat1
(lat1);
199
G_set_geodesic_distance_lat2
(lat2);
200
return
G_geodesic_distance_lon_to_lon
(lon1, lon2);
201
}
pi.h
G_geodesic_distance_lon_to_lon
double G_geodesic_distance_lon_to_lon(double lon1, double lon2)
Calculates geodesic distance.
Definition:
geodist.c:121
G_begin_geodesic_distance
void G_begin_geodesic_distance(double a, double e2)
Begin geodesic distance.
Definition:
geodist.c:50
G_set_geodesic_distance_lat2
void G_set_geodesic_distance_lat2(double lat2)
Sets geodesic distance lat2.
Definition:
geodist.c:83
x
#define x
t
double t
Definition:
driver/set_window.c:5
state
struct state state
Definition:
parser.c:101
G_set_geodesic_distance_lat1
void G_set_geodesic_distance_lat1(double lat1)
Sets geodesic distance lat1.
Definition:
geodist.c:69
G_geodesic_distance
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition:
geodist.c:196
Radians
#define Radians(x)
Definition:
pi.h:6
st
struct state * st
Definition:
parser.c:102
gis
geodist.c
Generated on Tue Mar 24 2020 14:13:52 for GRASS GIS 7 Programmer's Manual by
1.8.17