GRASS GIS 8 Programmer's Manual  8.5.0dev(2025)-ebe6f0621d
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 void G_begin_geodesic_distance(double a, double e2)
50 {
51  st->al = a;
52  st->boa = sqrt(1 - e2);
53  st->f = 1 - st->boa;
54  st->ff64 = st->f * st->f / 64;
55 }
56 
57 /*!
58  * \brief Sets geodesic distance lat1.
59  *
60  * Set the first latitude.
61  *
62  * <b>Note:</b> Must be called first.
63  *
64  * \param[in] lat1 first latitude
65  * \return
66  */
68 {
69  st->t1r = atan(st->boa * tan(Radians(lat1)));
70 }
71 
72 /*!
73  * \brief Sets geodesic distance lat2.
74  *
75  * Set the second latitude.
76  *
77  * <b>Note:</b> Must be called second.
78  *
79  * \param[in] lat2 second latitude
80  */
82 {
83  double stm, ctm, sdtm, cdtm;
84  double tm, dtm;
85 
86  st->t2r = atan(st->boa * tan(Radians(lat2)));
87 
88  tm = (st->t1r + st->t2r) / 2;
89  dtm = (st->t2r - st->t1r) / 2;
90 
91  stm = sin(tm);
92  ctm = cos(tm);
93  sdtm = sin(dtm);
94  cdtm = cos(dtm);
95 
96  st->t1 = stm * cdtm;
97  st->t1 = st->t1 * st->t1 * 2;
98 
99  st->t2 = sdtm * ctm;
100  st->t2 = st->t2 * st->t2 * 2;
101 
102  st->t3 = sdtm * sdtm;
103  st->t4 = cdtm * cdtm - stm * stm;
104 }
105 
106 /*!
107  * \brief Calculates geodesic distance.
108  *
109  * Calculates the geodesic distance from <i>lon1,lat1</i> to
110  * <i>lon2,lat2</i> in meters where <i>lat1</i> was the latitude
111  * passed to G_set_geodesic_distance_latl() and <i>lat2</i> was the
112  * latitude passed to G_set_geodesic_distance_lat2().
113  *
114  * \param[in] lon1 first longitude
115  * \param[in] lon2 second longitude
116  *
117  * \return double distance in meters
118  */
119 double G_geodesic_distance_lon_to_lon(double lon1, double lon2)
120 {
121  double a, cd, d, e, /*dl, */
122  q, sd, sdlmr, t, u, v, x, y;
123 
124  sdlmr = sin(Radians(lon2 - lon1) / 2);
125 
126  /* special case - shapiro */
127  if (sdlmr == 0.0 && st->t1r == st->t2r)
128  return 0.0;
129 
130  q = st->t3 + sdlmr * sdlmr * st->t4;
131 
132  /* special case - shapiro */
133  if (q == 1.0)
134  return M_PI * st->al;
135 
136  /* Mod: shapiro
137  * cd=1-2q is ill-conditioned if q is small O(10**-23)
138  * (for high lats? with lon1-lon2 < .25 degrees?)
139  * the computation of cd = 1-2*q will give cd==1.0.
140  * However, note that t=dl/sd is dl/sin(dl) which approaches 1 as dl->0.
141  * So the first step is to compute a good value for sd without using sin()
142  * and then check cd && q to see if we got cd==1.0 when we shouldn't.
143  * Note that dl isn't used except to get t,
144  * but both cd and sd are used later
145  */
146 
147  /* original code
148  cd=1-2*q;
149  dl=acos(cd);
150  sd=sin(dl);
151  t=dl/sd;
152  */
153 
154  cd = 1 - 2 * q; /* ill-conditioned subtraction for small q */
155  /* mod starts here */
156  sd = 2 * sqrt(q - q * q); /* sd^2 = 1 - cd^2 */
157  if (q != 0.0 && cd == 1.0) /* test for small q */
158  t = 1.0;
159  else if (sd == 0.0)
160  t = 1.0;
161  else
162  t = acos(cd) / sd; /* don't know how to fix acos(1-2*q) yet */
163  /* mod ends here */
164 
165  u = st->t1 / (1 - q);
166  v = st->t2 / q;
167  d = 4 * t * t;
168  x = u + v;
169  e = -2 * cd;
170  y = u - v;
171  a = -d * e;
172 
173  return st->al * sd *
174  (t - st->f / 4 * (t * x - y) +
175  st->ff64 * (x * (a + (t - (a + e) / 2) * x) + y * (-2 * d + e * y) +
176  d * x * y));
177 }
178 
179 /*!
180  * \brief Calculates geodesic distance.
181  *
182  * Calculates the geodesic distance from <i>lon1,lat1</i> to
183  * <i>lon2,lat2</i> in meters.
184  *
185  * <b>Note:</b> The calculation of the geodesic distance is fairly
186  * costly.
187  *
188  * \param[in] lon1,lat1 longitude,latitude of first point
189  * \param[in] lon2,lat2 longitude,latitude of second point
190  *
191  * \return distance in meters
192  */
193 double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
194 {
197  return G_geodesic_distance_lon_to_lon(lon1, lon2);
198 }
void G_set_geodesic_distance_lat2(double lat2)
Sets geodesic distance lat2.
Definition: geodist.c:81
void G_set_geodesic_distance_lat1(double lat1)
Sets geodesic distance lat1.
Definition: geodist.c:67
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition: geodist.c:193
double G_geodesic_distance_lon_to_lon(double lon1, double lon2)
Calculates geodesic distance.
Definition: geodist.c:119
void G_begin_geodesic_distance(double a, double e2)
Begin geodesic distance.
Definition: geodist.c:49
#define M_PI
Definition: gis.h:157
struct state state
Definition: parser.c:103
struct state * st
Definition: parser.c:104
#define Radians(x)
Definition: pi.h:6
double t
Definition: r_raster.c:39
#define x