geodesic.h File Reference

Header for the geodesic routines in C. More...

Go to the source code of this file.

Classes
struct	geod_geodesic
struct	geod_geodesicline
struct	geod_polygon

Macros
#define	GEODESIC_VERSION_MAJOR   1
#define	GEODESIC_VERSION_MINOR   32
#define	GEODESIC_VERSION_PATCH   0

Enumerations
enum	geod_mask {   GEOD_NONE = 0U , GEOD_LATITUDE = 1U << 7 | 0U , GEOD_LONGITUDE = 1U << 8 | 1U << 3 , GEOD_AZIMUTH = 1U << 9 | 0U ,   GEOD_DISTANCE = 1U << 10 | 1U << 0 , GEOD_DISTANCE_IN = 1U << 11 | 1U << 0 | 1U << 1 , GEOD_REDUCEDLENGTH = 1U << 12 | 1U << 0 | 1U << 2 , GEOD_GEODESICSCALE = 1U << 13 | 1U << 0 | 1U << 2 ,   GEOD_AREA = 1U << 14 | 1U << 4 , GEOD_ALL = 0x7F80U | 0x1FU }

Functions
void	geod_init (struct geod_geodesic *g, double a, double f)
void	geod_lineinit (struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)
void	geod_direct (const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, double *plat2, double *plon2, double *pazi2)
void	geod_inverse (const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)
void	geod_position (const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)
double	geod_gendirect (const struct geod_geodesic *g, double lat1, double lon1, double azi1, int arcmode, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)
double	geod_geninverse (const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2, double *pm12, double *pM12, double *pM21, double *pS12)
double	geod_genposition (const struct geod_geodesicline *l, int arcmode, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)
void	geod_polygon_init (struct geod_polygon *p, int polylinep)
void	geod_polygon_addpoint (const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)
void	geod_polygon_addedge (const struct geod_geodesic *g, struct geod_polygon *p, double azi, double s)
unsigned	geod_polygon_compute (const struct geod_geodesic *g, const struct geod_polygon *p, int reverse, int sign, double *pA, double *pP)
unsigned	geod_polygon_testpoint (const struct geod_geodesic *g, const struct geod_polygon *p, double lat, double lon, int reverse, int sign, double *pA, double *pP)
unsigned	geod_polygon_testedge (const struct geod_geodesic *g, const struct geod_polygon *p, double azi, double s, int reverse, int sign, double *pA, double *pP)
void	geod_polygonarea (const struct geod_geodesic *g, double lats[], double lons[], int n, double *pA, double *pP)

Detailed Description

Header for the geodesic routines in C.

This an implementation in C of the geodesic algorithms described in

• C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43–55 (2013); DOI: 10.1007/s00190-012-0578-z; addenda: geod-addenda.html.

The principal advantages of these algorithms over previous ones (e.g., Vincenty, 1975) are

• accurate to round off for |f| < 1/50;
• the solution of the inverse problem is always found;
• differential and integral properties of geodesics are computed.

The shortest path between two points on the ellipsoid at (lat1, lon1) and (lat2, lon2) is called the geodesic. Its length is s12 and the geodesic from point 1 to point 2 has forward azimuths azi1 and azi2 at the two end points.

Traditionally two geodesic problems are considered:

• the direct problem – given lat1, lon1, s12, and azi1, determine lat2, lon2, and azi2. This is solved by the function geod_direct().
• the inverse problem – given lat1, lon1, and lat2, lon2, determine s12, azi1, and azi2. This is solved by the function geod_inverse().

The ellipsoid is specified by its equatorial radius a (typically in meters) and flattening f. The routines are accurate to round off with double precision arithmetic provided that |f| < 1/50; for the WGS84 ellipsoid, the errors are less than 15 nanometers. (Reasonably accurate results are obtained for |f| < 1/5.) For a prolate ellipsoid, specify f < 0.

The routines also calculate several other quantities of interest

• S12 is the area between the geodesic from point 1 to point 2 and the equator; i.e., it is the area, measured counter-clockwise, of the quadrilateral with corners (lat1,lon1), (0,lon1), (0,lon2), and (lat2,lon2).
• m12, the reduced length of the geodesic is defined such that if the initial azimuth is perturbed by dazi1 (radians) then the second point is displaced by m12 dazi1 in the direction perpendicular to the geodesic. On a curved surface the reduced length obeys a symmetry relation, m12 + m21 = 0. On a flat surface, we have m12 = s12.
• M12 and M21 are geodesic scales. If two geodesics are parallel at point 1 and separated by a small distance dt, then they are separated by a distance M12 dt at point 2. M21 is defined similarly (with the geodesics being parallel to one another at point 2). On a flat surface, we have M12 = M21 = 1.
• a12 is the arc length on the auxiliary sphere. This is a construct for converting the problem to one in spherical trigonometry. a12 is measured in degrees. The spherical arc length from one equator crossing to the next is always 180°.

If points 1, 2, and 3 lie on a single geodesic, then the following addition rules hold:

• s13 = s12 + s23
• a13 = a12 + a23
• S13 = S12 + S23
• m13 = m12 M23 + m23 M21
• M13 = M12 M23 − (1 − M12 M21) m23 / m12
• M31 = M32 M21 − (1 − M23 M32) m12 / m23

The shortest distance returned by the solution of the inverse problem is (obviously) uniquely defined. However, in a few special cases there are multiple azimuths which yield the same shortest distance. Here is a catalog of those cases:

• lat1 = −lat2 (with neither at a pole). If azi1 = azi2, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting [azi1, azi2] = [azi2, azi1], [M12, M21] = [M21, M12], S12 = −S12. (This occurs when the longitude difference is near ±180° for oblate ellipsoids.)
• lon2 = lon1 ± 180° (with neither at a pole). If azi1 = 0° or ±180°, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting [azi1, azi2] = [−azi1, −azi2], S12 = −S12. (This occurs when the lat2 is near −lat1 for prolate ellipsoids.)
• Points 1 and 2 at opposite poles. There are infinitely many geodesics which can be generated by setting [azi1, azi2] = [azi1, azi2] + [d, −d], for arbitrary d. (For spheres, this prescription applies when points 1 and 2 are antipodal.)
• s12 = 0 (coincident points). There are infinitely many geodesics which can be generated by setting [azi1, azi2] = [azi1, azi2] + [d, d], for arbitrary d.

These routines are a simple transcription of the corresponding C++ classes in GeographicLib. The "class data" is represented by the structs geod_geodesic, geod_geodesicline, geod_polygon and pointers to these objects are passed as initial arguments to the member functions. Most of the internal comments have been retained. However, in the process of transcription some documentation has been lost and the documentation for the C++ classes, GeographicLib::Geodesic, GeographicLib::GeodesicLine, and GeographicLib::PolygonArea, should be consulted. The C++ code remains the "reference implementation". Think twice about restructuring the internals of the C code since this may make porting fixes from the C++ code more difficult.

Copyright (c) Charles Karney (2012-2013) charl.nosp@m.es@k.nosp@m.arney.nosp@m..com and licensed under the MIT/X11 License. For more information, see http://geographiclib.sourceforge.net/

This library was distributed with GeographicLib 1.32.

Macro Definition Documentation

GEODESIC_VERSION_MAJOR

#define GEODESIC_VERSION_MAJOR   1

The major version of the geodesic library. (This tracks the version of GeographicLib.)

GEODESIC_VERSION_MINOR

#define GEODESIC_VERSION_MINOR   32

The minor version of the geodesic library. (This tracks the version of GeographicLib.)

GEODESIC_VERSION_PATCH

#define GEODESIC_VERSION_PATCH   0

The patch level of the geodesic library. (This tracks the version of GeographicLib.)

Enumeration Type Documentation

geod_mask

enum geod_mask

mask values for the the caps argument to geod_lineinit().

Enumerator
GEOD_NONE	Calculate nothing
GEOD_LATITUDE	Calculate latitude
GEOD_LONGITUDE	Calculate longitude
GEOD_AZIMUTH	Calculate azimuth
GEOD_DISTANCE	Calculate distance
GEOD_DISTANCE_IN	Allow distance as input
GEOD_REDUCEDLENGTH	Calculate reduced length
GEOD_GEODESICSCALE	Calculate geodesic scale
GEOD_AREA	Calculate reduced length
GEOD_ALL	Calculate everything

Function Documentation

geod_direct()

void geod_direct	(	const struct geod_geodesic *	g,
		double	lat1,
		double	lon1,
		double	azi1,
		double	s12,
		double *	plat2,
		double *	plon2,
		double *	pazi2
	)

Solve the direct geodesic problem.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in]	lat1	latitude of point 1 (degrees).
[in]	lon1	longitude of point 1 (degrees).
[in]	azi1	azimuth at point 1 (degrees).
[in]	s12	distance between point 1 and point 2 (meters); it can be negative.
[out]	plat2	pointer to the latitude of point 2 (degrees).
[out]	plon2	pointer to the longitude of point 2 (degrees).
[out]	pazi2	pointer to the (forward) azimuth at point 2 (degrees).

g must have been initialized with a call to geod_init(). lat1 should be in the range [−90°, 90°]; lon1 and azi1 should be in the range [−540°, 540°). The values of lon2 and azi2 returned are in the range [−180°, 180°). Any of the "return" arguments plat2, etc., may be replaced by 0, if you do not need some quantities computed.

If either point is at a pole, the azimuth is defined by keeping the longitude fixed, writing lat = ±(90° − ε), and taking the limit ε → 0+. An arc length greater that 180° signifies a geodesic which is not a shortest path. (For a prolate ellipsoid, an additional condition is necessary for a shortest path: the longitudinal extent must not exceed of 180°.)

Example, determine the point 10000 km NE of JFK:

struct geod_geodesic g;
double lat, lon;
geod_init(&g, 6378137, 1/298.257223563);
geod_direct(&g, 40.64, -73.78, 45.0, 10e6, &lat, &lon, 0);
printf("%.5f %.5f\n", lat, lon);

geod_gendirect()

double geod_gendirect	(	const struct geod_geodesic *	g,
		double	lat1,
		double	lon1,
		double	azi1,
		int	arcmode,
		double	s12_a12,
		double *	plat2,
		double *	plon2,
		double *	pazi2,
		double *	ps12,
		double *	pm12,
		double *	pM12,
		double *	pM21,
		double *	pS12
	)

The general direct geodesic problem.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in]	lat1	latitude of point 1 (degrees).
[in]	lon1	longitude of point 1 (degrees).
[in]	azi1	azimuth at point 1 (degrees).
[in]	arcmode	flag determining the meaning of the s12_a12.
[in]	s12_a12	if arcmode is 0, this is the distance between point 1 and point 2 (meters); otherwise it is the arc length between point 1 and point 2 (degrees); it can be negative.
[out]	plat2	pointer to the latitude of point 2 (degrees).
[out]	plon2	pointer to the longitude of point 2 (degrees).
[out]	pazi2	pointer to the (forward) azimuth at point 2 (degrees).
[out]	ps12	pointer to the distance between point 1 and point 2 (meters).
[out]	pm12	pointer to the reduced length of geodesic (meters).
[out]	pM12	pointer to the geodesic scale of point 2 relative to point 1 (dimensionless).
[out]	pM21	pointer to the geodesic scale of point 1 relative to point 2 (dimensionless).
[out]	pS12	pointer to the area under the geodesic (meters2).

Returns

a12 arc length of between point 1 and point 2 (degrees).

g must have been initialized with a call to geod_init(). lat1 should be in the range [−90°, 90°]; lon1 and azi1 should be in the range [−540°, 540°). The function value a12 equals s12_a12 is arcmode is non-zero. Any of the "return" arguments plat2, etc., may be replaced by 0, if you do not need some quantities computed.

geod_geninverse()

double geod_geninverse	(	const struct geod_geodesic *	g,
		double	lat1,
		double	lon1,
		double	lat2,
		double	lon2,
		double *	ps12,
		double *	pazi1,
		double *	pazi2,
		double *	pm12,
		double *	pM12,
		double *	pM21,
		double *	pS12
	)

The general inverse geodesic calculation.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in]	lat1	latitude of point 1 (degrees).
[in]	lon1	longitude of point 1 (degrees).
[in]	lat2	latitude of point 2 (degrees).
[in]	lon2	longitude of point 2 (degrees).
[out]	ps12	pointer to the distance between point 1 and point 2 (meters).
[out]	pazi1	pointer to the azimuth at point 1 (degrees).
[out]	pazi2	pointer to the (forward) azimuth at point 2 (degrees).
[out]	pm12	pointer to the reduced length of geodesic (meters).
[out]	pM12	pointer to the geodesic scale of point 2 relative to point 1 (dimensionless).
[out]	pM21	pointer to the geodesic scale of point 1 relative to point 2 (dimensionless).
[out]	pS12	pointer to the area under the geodesic (meters2).

Returns

a12 arc length of between point 1 and point 2 (degrees).

g must have been initialized with a call to geod_init(). lat1 and lat2 should be in the range [−90°, 90°]; lon1 and lon2 should be in the range [−540°, 540°). Any of the "return" arguments ps12, etc., may be replaced by 0, if you do not need some quantities computed.

geod_genposition()

double geod_genposition	(	const struct geod_geodesicline *	l,
		int	arcmode,
		double	s12_a12,
		double *	plat2,
		double *	plon2,
		double *	pazi2,
		double *	ps12,
		double *	pm12,
		double *	pM12,
		double *	pM21,
		double *	pS12
	)

The general position function.

Parameters

[in]	l	a pointer to the geod_geodesicline object specifying the geodesic line.
[in]	arcmode	flag determining the meaning of the second parameter; if arcmode is 0, then l must have been initialized with caps |= GEOD_DISTANCE_IN.
[in]	s12_a12	if arcmode is 0, this is the distance between point 1 and point 2 (meters); otherwise it is the arc length between point 1 and point 2 (degrees); it can be negative.
[out]	plat2	pointer to the latitude of point 2 (degrees).
[out]	plon2	pointer to the longitude of point 2 (degrees); requires that l was initialized with caps |= GEOD_LONGITUDE.
[out]	pazi2	pointer to the (forward) azimuth at point 2 (degrees).
[out]	ps12	pointer to the distance between point 1 and point 2 (meters); requires that l was initialized with caps |= GEOD_DISTANCE.
[out]	pm12	pointer to the reduced length of geodesic (meters); requires that l was initialized with caps |= GEOD_REDUCEDLENGTH.
[out]	pM12	pointer to the geodesic scale of point 2 relative to point 1 (dimensionless); requires that l was initialized with caps |= GEOD_GEODESICSCALE.
[out]	pM21	pointer to the geodesic scale of point 1 relative to point 2 (dimensionless); requires that l was initialized with caps |= GEOD_GEODESICSCALE.
[out]	pS12	pointer to the area under the geodesic (meters2); requires that l was initialized with caps |= GEOD_AREA.

Returns

a12 arc length of between point 1 and point 2 (degrees).

l must have been initialized with a call to geod_lineinit() with caps |= GEOD_DISTANCE_IN. The values of lon2 and azi2 returned are in the range [−180°, 180°). Any of the "return" arguments plat2, etc., may be replaced by 0, if you do not need some quantities computed. Requesting a value which l is not capable of computing is not an error; the corresponding argument will not be altered.

Example, compute way points between JFK and Singapore Changi Airport using geod_genposition(). In this example, the points are evenly space in arc length (and so only approximately equally space in distance). This is faster than using geod_position() would be appropriate if drawing the path on a map.

struct geod_geodesic g;
struct geod_geodesicline l;
double a12, azi1, lat[101],lon[101];
int i;
geod_init(&g, 6378137, 1/298.257223563);
a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99,
                      0, &azi1, 0, 0, 0, 0, 0);
geod_lineinit(&l, &g, 40.64, -73.78, azi1, GEOD_LATITUDE | GEOD_LONGITUDE);
for (i = 0; i < 101; ++i) {
  geod_genposition(&l, 1, i * a12 * 0.01,
                   lat + i, lon + i, 0, 0, 0, 0, 0, 0);
  printf("%.5f %.5f\n", lat[i], lon[i]);
}

geod_init()

void geod_init	(	struct geod_geodesic *	g,
		double	a,
		double	f
	)

Initialize a geod_geodesic object.

Parameters

[out]	g	a pointer to the object to be initialized.
[in]	a	the equatorial radius (meters).
[in]	f	the flattening.

geod_inverse()

void geod_inverse	(	const struct geod_geodesic *	g,
		double	lat1,
		double	lon1,
		double	lat2,
		double	lon2,
		double *	ps12,
		double *	pazi1,
		double *	pazi2
	)

Solve the inverse geodesic problem.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in]	lat1	latitude of point 1 (degrees).
[in]	lon1	longitude of point 1 (degrees).
[in]	lat2	latitude of point 2 (degrees).
[in]	lon2	longitude of point 2 (degrees).
[out]	ps12	pointer to the distance between point 1 and point 2 (meters).
[out]	pazi1	pointer to the azimuth at point 1 (degrees).
[out]	pazi2	pointer to the (forward) azimuth at point 2 (degrees).

g must have been initialized with a call to geod_init(). lat1 and lat2 should be in the range [−90°, 90°]; lon1 and lon2 should be in the range [−540°, 540°). The values of azi1 and azi2 returned are in the range [−180°, 180°). Any of the "return" arguments ps12, etc., may be replaced by 0, if you do not need some quantities computed.

If either point is at a pole, the azimuth is defined by keeping the longitude fixed, writing lat = ±(90° − ε), and taking the limit ε → 0+.

The solution to the inverse problem is found using Newton's method. If this fails to converge (this is very unlikely in geodetic applications but does occur for very eccentric ellipsoids), then the bisection method is used to refine the solution.

Example, determine the distance between JFK and Singapore Changi Airport:

struct geod_geodesic g;
double s12;
geod_init(&g, 6378137, 1/298.257223563);
geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, 0, 0);
printf("%.3f\n", s12);

geod_lineinit()

void geod_lineinit	(	struct geod_geodesicline *	l,
		const struct geod_geodesic *	g,
		double	lat1,
		double	lon1,
		double	azi1,
		unsigned	caps
	)

Initialize a geod_geodesicline object.

Parameters

[out]	l	a pointer to the object to be initialized.
[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in]	lat1	latitude of point 1 (degrees).
[in]	lon1	longitude of point 1 (degrees).
[in]	azi1	azimuth at point 1 (degrees).
[in]	caps	bitor'ed combination of geod_mask() values specifying the capabilities the geod_geodesicline object should possess, i.e., which quantities can be returned in calls to geod_position() and geod_genposition().

g must have been initialized with a call to geod_init(). lat1 should be in the range [−90°, 90°]; lon1 and azi1 should be in the range [−540°, 540°).

The geod_mask values are [see geod_mask()]:

• caps |= GEOD_LATITUDE for the latitude lat2; this is added automatically,
• caps |= GEOD_LONGITUDE for the latitude lon2,
• caps |= GEOD_AZIMUTH for the latitude azi2; this is added automatically,
• caps |= GEOD_DISTANCE for the distance s12,
• caps |= GEOD_REDUCEDLENGTH for the reduced length m12,
• caps |= GEOD_GEODESICSCALE for the geodesic scales M12 and M21,
• caps |= GEOD_AREA for the area S12,
• caps |= GEOD_DISTANCE_IN permits the length of the geodesic to be given in terms of s12; without this capability the length can only be specified in terms of arc length.

A value of caps = 0 is treated as GEOD_LATITUDE | GEOD_LONGITUDE | GEOD_AZIMUTH | GEOD_DISTANCE_IN (to support the solution of the "standard" direct problem).

geod_polygon_addedge()

void geod_polygon_addedge	(	const struct geod_geodesic *	g,
		struct geod_polygon *	p,
		double	azi,
		double	s
	)

Add an edge to the polygon or polyline.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in,out]	p	a pointer to the geod_polygon object specifying the polygon.
[in]	azi	azimuth at current point (degrees).
[in]	s	distance from current point to next point (meters).

g and p must have been initialized with calls to geod_init() and geod_polygon_init(), respectively. The same g must be used for all the points and edges in a polygon. azi should be in the range [−540°, 540°). This does nothing if no points have been added yet. The lat and lon fields of p give the location of the new vertex.

geod_polygon_addpoint()

void geod_polygon_addpoint	(	const struct geod_geodesic *	g,
		struct geod_polygon *	p,
		double	lat,
		double	lon
	)

Add a point to the polygon or polyline.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in,out]	p	a pointer to the geod_polygon object specifying the polygon.
[in]	lat	the latitude of the point (degrees).
[in]	lon	the longitude of the point (degrees).

g and p must have been initialized with calls to geod_init() and geod_polygon_init(), respectively. The same g must be used for all the points and edges in a polygon. lat should be in the range [−90°, 90°] and lon should be in the range [−540°, 540°).

An example of the use of this function is given in the documentation for geod_polygon_compute().

geod_polygon_compute()

unsigned geod_polygon_compute	(	const struct geod_geodesic *	g,
		const struct geod_polygon *	p,
		int	reverse,
		int	sign,
		double *	pA,
		double *	pP
	)

Return the results for a polygon.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in]	p	a pointer to the geod_polygon object specifying the polygon.
[in]	reverse	if non-zero then clockwise (instead of counter-clockwise) traversal counts as a positive area.
[in]	sign	if non-zero then return a signed result for the area if the polygon is traversed in the "wrong" direction instead of returning the area for the rest of the earth.
[out]	pA	pointer to the area of the polygon (meters2); only set if polyline is non-zero in the call to geod_polygon_init().
[out]	pP	pointer to the perimeter of the polygon or length of the polyline (meters).

Returns

the number of points.

Only simple polygons (which are not self-intersecting) are allowed. There's no need to "close" the polygon by repeating the first vertex. Set pA or pP to zero, if you do not want the corresponding quantity returned.

Example, compute the perimeter and area of the geodesic triangle with vertices (0°N,0°E), (0°N,90°E), (90°N,0°E).

double A, P;
int n;
struct geod_geodesic g;
struct geod_polygon p;
geod_init(&g, 6378137, 1/298.257223563);
geod_polygon_init(&p, 0);
 
geod_polygon_addpoint(&g, &p,  0,  0);
geod_polygon_addpoint(&g, &p,  0, 90);
geod_polygon_addpoint(&g, &p, 90,  0);
n = geod_polygon_compute(&g, &p, 0, 1, &A, &P);
printf("%d %.8f %.3f\n", n, P, A);

geod_polygon_init()

void geod_polygon_init	(	struct geod_polygon *	p,
		int	polylinep
	)

Initialize a geod_polygon object.

Parameters

[out]	p	a pointer to the object to be initialized.
[in]	polylinep	non-zero if a polyline instead of a polygon.

If polylinep is zero, then the sequence of vertices and edges added by geod_polygon_addpoint() and geod_polygon_addedge() define a polygon and the perimeter and area are returned by geod_polygon_compute(). If polylinep is non-zero, then the vertices and edges define a polyline and only the perimeter is returned by geod_polygon_compute().

An example of the use of this function is given in the documentation for geod_polygon_compute().

geod_polygon_testedge()

unsigned geod_polygon_testedge	(	const struct geod_geodesic *	g,
		const struct geod_polygon *	p,
		double	azi,
		double	s,
		int	reverse,
		int	sign,
		double *	pA,
		double *	pP
	)

Return the results assuming a tentative final test point is added via an azimuth and distance; however, the data for the test point is not saved. This lets you report a running result for the perimeter and area as the user moves the mouse cursor. Ordinary floating point arithmetic is used to accumulate the data for the test point; thus the area and perimeter returned are less accurate than if geod_polygon_addedge() and geod_polygon_compute() are used.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in]	p	a pointer to the geod_polygon object specifying the polygon.
[in]	azi	azimuth at current point (degrees).
[in]	s	distance from current point to final test point (meters).
[in]	reverse	if non-zero then clockwise (instead of counter-clockwise) traversal counts as a positive area.
[in]	sign	if non-zero then return a signed result for the area if the polygon is traversed in the "wrong" direction instead of returning the area for the rest of the earth.
[out]	pA	pointer to the area of the polygon (meters2); only set if polyline is non-zero in the call to geod_polygon_init().
[out]	pP	pointer to the perimeter of the polygon or length of the polyline (meters).

Returns

the number of points.

azi should be in the range [−540°, 540°).

geod_polygon_testpoint()

unsigned geod_polygon_testpoint	(	const struct geod_geodesic *	g,
		const struct geod_polygon *	p,
		double	lat,
		double	lon,
		int	reverse,
		int	sign,
		double *	pA,
		double *	pP
	)

Return the results assuming a tentative final test point is added; however, the data for the test point is not saved. This lets you report a running result for the perimeter and area as the user moves the mouse cursor. Ordinary floating point arithmetic is used to accumulate the data for the test point; thus the area and perimeter returned are less accurate than if geod_polygon_addpoint() and geod_polygon_compute() are used.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in]	p	a pointer to the geod_polygon object specifying the polygon.
[in]	lat	the latitude of the test point (degrees).
[in]	lon	the longitude of the test point (degrees).
[in]	reverse	if non-zero then clockwise (instead of counter-clockwise) traversal counts as a positive area.
[in]	sign	if non-zero then return a signed result for the area if the polygon is traversed in the "wrong" direction instead of returning the area for the rest of the earth.
[out]	pA	pointer to the area of the polygon (meters2); only set if polyline is non-zero in the call to geod_polygon_init().
[out]	pP	pointer to the perimeter of the polygon or length of the polyline (meters).

Returns

the number of points.

lat should be in the range [−90°, 90°] and lon should be in the range [−540°, 540°).

geod_polygonarea()

void geod_polygonarea	(	const struct geod_geodesic *	g,
		double	lats[],
		double	lons[],
		int	n,
		double *	pA,
		double *	pP
	)

A simple interface for computing the area of a geodesic polygon.

Parameters

[in]	g	a pointer to the geod_geodesic object specifying the ellipsoid.
[in]	lats	an array of latitudes of the polygon vertices (degrees).
[in]	lons	an array of longitudes of the polygon vertices (degrees).
[in]	n	the number of vertices.
[out]	pA	pointer to the area of the polygon (meters2).
[out]	pP	pointer to the perimeter of the polygon (meters).

lats should be in the range [−90°, 90°]; lons should be in the range [−540°, 540°).

Only simple polygons (which are not self-intersecting) are allowed. There's no need to "close" the polygon by repeating the first vertex. The area returned is signed with counter-clockwise traversal being treated as positive.

Example, compute the area of Antarctic:

double
  lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
            -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7},
  lons[] = {-74, -102, -102, -131, -163, 163, 172, 140, 113,
             88, 59, 25, -4, -14, -33, -46, -61};
struct geod_geodesic g;
double A, P;
geod_init(&g, 6378137, 1/298.257223563);
geod_polygonarea(&g, lats, lons, (sizeof lats) / (sizeof lats[0]), &A, &P);
printf("%.0f %.2f\n", A, P);

geod_position()

void geod_position	(	const struct geod_geodesicline *	l,
		double	s12,
		double *	plat2,
		double *	plon2,
		double *	pazi2
	)

Compute the position along a geod_geodesicline.

Parameters

[in]	l	a pointer to the geod_geodesicline object specifying the geodesic line.
[in]	s12	distance between point 1 and point 2 (meters); it can be negative.
[out]	plat2	pointer to the latitude of point 2 (degrees).
[out]	plon2	pointer to the longitude of point 2 (degrees); requires that l was initialized with caps |= GEOD_LONGITUDE.
[out]	pazi2	pointer to the (forward) azimuth at point 2 (degrees).

l must have been initialized with a call to geod_lineinit() with caps |= GEOD_DISTANCE_IN. The values of lon2 and azi2 returned are in the range [−180°, 180°). Any of the "return" arguments plat2, etc., may be replaced by 0, if you do not need some quantities computed.

Example, compute way points between JFK and Singapore Changi Airport the "obvious" way using geod_direct():

struct geod_geodesic g;
double s12, azi1, lat[101],lon[101];
int i;
geod_init(&g, 6378137, 1/298.257223563);
geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
for (i = 0; i < 101; ++i) {
  geod_direct(&g, 40.64, -73.78, azi1, i * s12 * 0.01, lat + i, lon + i, 0);
  printf("%.5f %.5f\n", lat[i], lon[i]);
}

A faster way using geod_position():

struct geod_geodesic g;
struct geod_geodesicline l;
double s12, azi1, lat[101],lon[101];
int i;
geod_init(&g, 6378137, 1/298.257223563);
geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
geod_lineinit(&l, &g, 40.64, -73.78, azi1, 0);
for (i = 0; i < 101; ++i) {
  geod_position(&l, i * s12 * 0.01, lat + i, lon + i, 0);
  printf("%.5f %.5f\n", lat[i], lon[i]);
}