GRASS GIS 8 Programmer's Manual 8.2.1(2023)-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
27static 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
35static 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
50void 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
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 */
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 */
121double 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 */
196double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
197{
200 return G_geodesic_distance_lon_to_lon(lon1, lon2);
201}
double t
void G_set_geodesic_distance_lat2(double lat2)
Sets geodesic distance lat2.
Definition: geodist.c:83
void G_set_geodesic_distance_lat1(double lat1)
Sets geodesic distance lat1.
Definition: geodist.c:69
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition: geodist.c:196
double G_geodesic_distance_lon_to_lon(double lon1, double lon2)
Calculates geodesic distance.
Definition: geodist.c:121
void G_begin_geodesic_distance(double a, double e2)
Begin geodesic distance.
Definition: geodist.c:50
struct state state
Definition: parser.c:103
struct state * st
Definition: parser.c:104
#define Radians(x)
Definition: pi.h:6
#define x