GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-0306f9333d
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
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
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  */
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  sdlmr = sin(Radians(lon2 - lon1) / 2);
127
128  /* special case - shapiro */
129  if (sdlmr == 0.0 && st->t1r == st->t2r)
130  return 0.0;
131
132  q = st->t3 + sdlmr * sdlmr * st->t4;
133
134  /* special case - shapiro */
135  if (q == 1.0)
136  return M_PI * st->al;
137
138  /* Mod: shapiro
139  * cd=1-2q is ill-conditioned if q is small O(10**-23)
140  * (for high lats? with lon1-lon2 < .25 degrees?)
141  * the computation of cd = 1-2*q will give cd==1.0.
142  * However, note that t=dl/sd is dl/sin(dl) which approaches 1 as dl->0.
143  * So the first step is to compute a good value for sd without using sin()
144  * and then check cd && q to see if we got cd==1.0 when we shouldn't.
145  * Note that dl isn't used except to get t,
146  * but both cd and sd are used later
147  */
148
149  /* original code
150  cd=1-2*q;
151  dl=acos(cd);
152  sd=sin(dl);
153  t=dl/sd;
154  */
155
156  cd = 1 - 2 * q; /* ill-conditioned subtraction for small q */
157  /* mod starts here */
158  sd = 2 * sqrt(q - q * q); /* sd^2 = 1 - cd^2 */
159  if (q != 0.0 && cd == 1.0) /* test for small q */
160  t = 1.0;
161  else if (sd == 0.0)
162  t = 1.0;
163  else
164  t = acos(cd) / sd; /* don't know how to fix acos(1-2*q) yet */
165  /* mod ends here */
166
167  u = st->t1 / (1 - q);
168  v = st->t2 / q;
169  d = 4 * t * t;
170  x = u + v;
171  e = -2 * cd;
172  y = u - v;
173  a = -d * e;
174
175  return st->al * sd *
176  (t - st->f / 4 * (t * x - y) +
177  st->ff64 * (x * (a + (t - (a + e) / 2) * x) + y * (-2 * d + e * y) +
178  d * x * y));
179 }
180
181 /*!
182  * \brief Calculates geodesic distance.
183  *
184  * Calculates the geodesic distance from <i>lon1,lat1</i> to
185  * <i>lon2,lat2</i> in meters.
186  *
187  * <b>Note:</b> The calculation of the geodesic distance is fairly
188  * costly.
189  *
190  * \param lon1,lat1 longitude,latitude of first point
191  * \param lon2,lat2 longitude,latitude of second point
192  *
193  * \return distance in meters
194  */
195 double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
196 {
199  return G_geodesic_distance_lon_to_lon(lon1, lon2);
200 }
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:195
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
#define M_PI
Definition: gis.h:158
struct state state
Definition: parser.c:103
struct state * st
Definition: parser.c:104