GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
gis/distance.c
Go to the documentation of this file.
1/*!
2 \file lib/gis/distance.c
3
4 \brief GIS Library - Distance calculation functions.
5
6 WARNING: this code is preliminary and may be changed,
7 including calling sequences to any of the functions
8 defined here.
9
10 (C) 2001-2009, 2011 by the GRASS Development Team
11
12 This program is free software under the GNU General Public License
13 (>=v2). Read the file COPYING that comes with GRASS for details.
14
15 \author Original author CERL
16 */
17
18#include <math.h>
19#include <grass/gis.h>
20#include <grass/glocale.h>
21
22static double min4(double, double, double, double);
23static double min2(double, double);
24
25static struct state {
26 int projection;
27 double factor;
28} state;
29
30static struct state *st = &state;
31
32/*!
33 \brief Begin distance calculations.
34
35 Initializes the distance calculations. It is used both for the
36 planimetric and latitude-longitude projections.
37
38 \return 0 if projection has no metrix (ie. imagery)
39 \return 1 if projection is planimetric
40 \return 2 if projection is latitude-longitude
41 */
43{
44 double a, e2;
45
46 st->factor = 1.0;
47 switch (st->projection = G_projection()) {
48 case PROJECTION_LL:
51 return 2;
52 default:
54 if (st->factor <= 0.0) {
55 st->factor = 1.0; /* assume meter grid */
56 return 0;
57 }
58 return 1;
59 }
60}
61
62/*!
63 \brief Returns distance in meters.
64
65 This routine computes the distance, in meters, from
66 <i>x1</i>,<i>y1</i> to <i>x2</i>,<i>y2</i>. If the projection is
67 latitude-longitude, this distance is measured along the
68 geodesic. Two routines perform geodesic distance calculations.
69
70 \param e1,n1 east-north coordinates of first point
71 \param e2,n2 east-north coordinates of second point
72
73 \return distance
74 */
75double G_distance(double e1, double n1, double e2, double n2)
76{
77 if (st->projection == PROJECTION_LL)
78 return G_geodesic_distance(e1, n1, e2, n2);
79 else
80 return st->factor * hypot(e1 - e2, n1 - n2);
81}
82
83/*!
84 \brief Returns distance between two line segments in meters.
85
86 \param ax1,ay1,ax2,ay2 first segment
87 \param bx1,by1,bx2,by2 second segment
88
89 \return distance value
90 */
91double G_distance_between_line_segments(double ax1, double ay1, double ax2,
92 double ay2, double bx1, double by1,
93 double bx2, double by2)
94{
95 double ra, rb;
96 double x, y;
97
98 /* if the segments intersect, then the distance is zero */
99 if (G_intersect_line_segments(ax1, ay1, ax2, ay2, bx1, by1, bx2, by2, &ra,
100 &rb, &x, &y) > 0)
101 return 0.0;
102 return min4(G_distance_point_to_line_segment(ax1, ay1, bx1, by1, bx2, by2),
103 G_distance_point_to_line_segment(ax2, ay2, bx1, by1, bx2, by2),
104 G_distance_point_to_line_segment(bx1, by1, ax1, ay1, ax2, ay2),
105 G_distance_point_to_line_segment(bx2, by2, ax1, ay1, ax2, ay2));
106}
107
108/*!
109 \brief Returns distance between a point and line segment in meters.
110
111 \param xp,yp point coordinates
112 \param x1,y1 segment point coordinates
113 \param x2,y2 segment point coordinates
114
115 \return distance
116 */
117double G_distance_point_to_line_segment(double xp, double yp, double x1,
118 double y1, double x2, double y2)
119{
120 double dx, dy;
121 double x, y;
122 double xq, yq, ra, rb;
123 int t;
124
125 /* define the perpendicular to the segment through the point */
126 dx = x1 - x2;
127 dy = y1 - y2;
128
129 if (dx == 0.0 && dy == 0.0)
130 return G_distance(x1, y1, xp, yp);
131
132 if (fabs(dy) > fabs(dx)) {
133 xq = xp + dy;
134 yq = (dx / dy) * (xp - xq) + yp;
135 }
136 else {
137 yq = yp + dx;
138 xq = (dy / dx) * (yp - yq) + xp;
139 }
140
141 /* find the intersection of the perpendicular with the segment */
142 t = G_intersect_line_segments(xp, yp, xq, yq, x1, y1, x2, y2, &ra, &rb, &x,
143 &y);
144 switch (t) {
145 case 0:
146 case 1:
147 break;
148 default:
149 /* parallel/colinear cases shouldn't occur with perpendicular lines */
150 G_warning(_("%s: shouldn't happen: "
151 "code=%d P=(%f,%f) S=(%f,%f)(%f,%f)"),
152 "G_distance_point_to_line_segment", t, xp, yp, x1, y1, x2,
153 y2);
154 return -1.0;
155 }
156
157 /* if x,y lies on the segment, then the distance is from to x,y */
158 if (rb >= 0 && rb <= 1.0)
159 return G_distance(x, y, xp, yp);
160
161 /* otherwise the distance is the short of the distances to the endpoints
162 * of the segment
163 */
164 return min2(G_distance(x1, y1, xp, yp), G_distance(x2, y2, xp, yp));
165}
166
167static double min4(double a, double b, double c, double d)
168{
169 return min2(min2(a, b), min2(c, d));
170}
171
172static double min2(double a, double b)
173{
174 return a < b ? a : b;
175}
double b
double t
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition geodist.c:195
void G_begin_geodesic_distance(double a, double e2)
Begin geodesic distance.
Definition geodist.c:50
int G_get_ellipsoid_parameters(double *a, double *e2)
get ellipsoid parameters
Definition get_ellipse.c:66
double G_distance_point_to_line_segment(double xp, double yp, double x1, double y1, double x2, double y2)
Returns distance between a point and line segment in meters.
int G_begin_distance_calculations(void)
Begin distance calculations.
double G_distance_between_line_segments(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2)
Returns distance between two line segments in meters.
double G_distance(double e1, double n1, double e2, double n2)
Returns distance in meters.
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
int G_intersect_line_segments(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2, double *ra, double *rb, double *x, double *y)
Definition intersect.c:84
int G_projection(void)
Query cartographic projection.
Definition proj1.c:32
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition proj3.c:146
#define x