Loading...
Searching...
No Matches
gis/distance.c
Go to the documentation of this file.
/*!
\file lib/gis/distance.c
\brief GIS Library - Distance calculation functions.
WARNING: this code is preliminary and may be changed,
including calling sequences to any of the functions
defined here.
(C) 2001-2009, 2011 by the GRASS Development Team
This program is free software under the GNU General Public License
(>=v2). Read the file COPYING that comes with GRASS for details.
\author Original author CERL
*/
#include <math.h>
#include <grass/gis.h>
#include <grass/glocale.h>
static double min4(double, double, double, double);
static double min2(double, double);
static struct state {
int projection;
double factor;
} state;
static struct state *st = &state;
/*!
\brief Begin distance calculations.
Initializes the distance calculations. It is used both for the
planimetric and latitude-longitude projections.
\return 0 if projection has no metrics (ie. imagery)
\return 1 if projection is planimetric
\return 2 if projection is latitude-longitude
*/
int G_begin_distance_calculations(void)
{
double a, e2;
st->factor = 1.0;
switch (st->projection = G_projection()) {
case PROJECTION_LL:
G_get_ellipsoid_parameters(&a, &e2);
G_begin_geodesic_distance(a, e2);
return 2;
default:
st->factor = G_database_units_to_meters_factor();
if (st->factor <= 0.0) {
st->factor = 1.0; /* assume meter grid */
return 0;
}
return 1;
}
}
/*!
\brief Returns distance in meters.
This routine computes the distance, in meters, from
<i>x1</i>,<i>y1</i> to <i>x2</i>,<i>y2</i>. If the projection is
latitude-longitude, this distance is measured along the
geodesic. Two routines perform geodesic distance calculations.
\param[in] e1,n1 east-north coordinates of first point
\param[in] e2,n2 east-north coordinates of second point
\return distance
*/
double G_distance(double e1, double n1, double e2, double n2)
{
if (st->projection == PROJECTION_LL)
return G_geodesic_distance(e1, n1, e2, n2);
else
return st->factor * hypot(e1 - e2, n1 - n2);
}
/*!
\brief Returns distance between two line segments in meters.
\param ax1,ay1,ax2,ay2 first segment
\param bx1,by1,bx2,by2 second segment
\return distance value
*/
double G_distance_between_line_segments(double ax1, double ay1, double ax2,
double ay2, double bx1, double by1,
double bx2, double by2)
{
double ra, rb;
double x, y;
/* if the segments intersect, then the distance is zero */
if (G_intersect_line_segments(ax1, ay1, ax2, ay2, bx1, by1, bx2, by2, &ra,
&rb, &x, &y) > 0)
return 0.0;
return min4(G_distance_point_to_line_segment(ax1, ay1, bx1, by1, bx2, by2),
G_distance_point_to_line_segment(ax2, ay2, bx1, by1, bx2, by2),
G_distance_point_to_line_segment(bx1, by1, ax1, ay1, ax2, ay2),
G_distance_point_to_line_segment(bx2, by2, ax1, ay1, ax2, ay2));
}
/*!
\brief Returns distance between a point and line segment in meters.
\param xp,yp point coordinates
\param x1,y1 segment point coordinates
\param x2,y2 segment point coordinates
\return distance
*/
double G_distance_point_to_line_segment(double xp, double yp, double x1,
double y1, double x2, double y2)
{
double dx, dy;
double x, y;
double xq, yq, ra, rb;
int t;
/* define the perpendicular to the segment through the point */
dx = x1 - x2;
dy = y1 - y2;
if (dx == 0.0 && dy == 0.0)
return G_distance(x1, y1, xp, yp);
if (fabs(dy) > fabs(dx)) {
xq = xp + dy;
yq = (dx / dy) * (xp - xq) + yp;
}
else {
yq = yp + dx;
xq = (dy / dx) * (yp - yq) + xp;
}
/* find the intersection of the perpendicular with the segment */
t = G_intersect_line_segments(xp, yp, xq, yq, x1, y1, x2, y2, &ra, &rb, &x,
&y);
switch (t) {
case 0:
case 1:
break;
default:
/* parallel/colinear cases shouldn't occur with perpendicular lines */
G_warning(_("%s: shouldn't happen: "
"code=%d P=(%f,%f) S=(%f,%f)(%f,%f)"),
"G_distance_point_to_line_segment", t, xp, yp, x1, y1, x2,
y2);
return -1.0;
}
/* if x,y lies on the segment, then the distance is from to x,y */
if (rb >= 0 && rb <= 1.0)
return G_distance(x, y, xp, yp);
/* otherwise the distance is the short of the distances to the endpoints
* of the segment
*/
return min2(G_distance(x1, y1, xp, yp), G_distance(x2, y2, xp, yp));
}
static double min4(double a, double b, double c, double d)
{
return min2(min2(a, b), min2(c, d));
}
static double min2(double a, double b)
{
return a < b ? a : b;
}
b
double b
Definition driver/set_window.c:5
t
double t
Definition driver/set_window.c:5
G_geodesic_distance
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition geodist.c:193
G_begin_geodesic_distance
void G_begin_geodesic_distance(double a, double e2)
Begin geodesic distance.
Definition geodist.c:49
G_get_ellipsoid_parameters
int G_get_ellipsoid_parameters(double *a, double *e2)
get ellipsoid parameters
Definition get_ellipse.c:66
G_distance_point_to_line_segment
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.
Definition gis/distance.c:117
G_begin_distance_calculations
int G_begin_distance_calculations(void)
Begin distance calculations.
Definition gis/distance.c:42
G_distance_between_line_segments
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.
Definition gis/distance.c:91
G_distance
double G_distance(double e1, double n1, double e2, double n2)
Returns distance in meters.
Definition gis/distance.c:75
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
G_intersect_line_segments
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
G_projection
int G_projection(void)
Query cartographic projection.
Definition proj1.c:32
G_database_units_to_meters_factor
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition proj3.c:146
x
#define x