GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71933
 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 
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 
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 {
200  return G_geodesic_distance_lon_to_lon(lon1, lon2);
201 }
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:115
#define x
struct state * st
Definition: parser.c:103
void G_set_geodesic_distance_lat2(double lat2)
Sets geodesic distance lat2.
Definition: geodist.c:83
double t
Definition: r_raster.c:39
#define Radians(x)
Definition: pi.h:6
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition: geodist.c:196
void G_set_geodesic_distance_lat1(double lat1)
Sets geodesic distance lat1.
Definition: geodist.c:69
struct state state
Definition: parser.c:102