C library for Geodesics: geodesic.h File Reference (original) (raw)

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

Go to the source code of this file.

API for the geodesic routines in C.

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::PolygonAreaT, 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-2022) charl.nosp@m.es@k.nosp@m.arney.nosp@m..com and licensed under the MIT/X11 License. For more information, see https://geographiclib.sourceforge.io/

Definition in file geodesic.h.

◆ GEODESIC_VERSION_MAJOR

#define GEODESIC_VERSION_MAJOR

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

Definition at line 29 of file geodesic.h.

◆ GEODESIC_VERSION_MINOR

#define GEODESIC_VERSION_MINOR

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

Definition at line 34 of file geodesic.h.

◆ GEODESIC_VERSION_PATCH

#define GEODESIC_VERSION_PATCH

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

Definition at line 39 of file geodesic.h.

◆ GEODESIC_VERSION_NUM

| #define GEODESIC_VERSION_NUM | ( | | a, | | ------------------------------ | - | | -- | | | b, | | | | | | c | | | | | ) | | | |

Pack the version components into a single integer. Users should not rely on this particular packing of the components of the version number; see the documentation for GEODESIC_VERSION, below.

Definition at line 46 of file geodesic.h.

◆ GEODESIC_VERSION

The version of the geodesic library as a single integer, packed as MMmmmmpp where MM is the major version, mmmm is the minor version, and pp is the patch level. Users should not rely on this particular packing of the components of the version number. Instead they should use a test such as

#if GEODESIC_VERSION >= GEODESIC_VERSION_NUM(1,40,0)

...

#endif

Definition at line 59 of file geodesic.h.

◆ GEOD_DLL

◆ geod_mask

mask values for the caps argument to geod_lineinit().

Definition at line 810 of file geodesic.h.

◆ geod_flags

flag values for the flags argument to geod_gendirect() and geod_genposition()

Definition at line 827 of file geodesic.h.

◆ geod_init()

Initialize a geod_geodesic object.

Parameters

Referenced by main().

◆ geod_direct()

Solve the direct geodesic problem.

Parameters

g must have been initialized with a call to geod_init(). lat1 should be in the range [−90°, 90°]. 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:

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);

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

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

Referenced by main().

◆ geod_gendirect()

The general direct geodesic problem.

Parameters

Returns

a12 arc length from point 1 to point 2 (degrees).

g must have been initialized with a call to geod_init(). lat1 should be in the range [−90°, 90°]. The function value a12 equals s12_a12 if flags & GEOD_ARCMODE. Any of the "return" arguments, plat2, etc., may be replaced by 0, if you do not need some quantities computed.

With flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so that the quantity lon2 − lon1 indicates how many times and in what sense the geodesic encircles the ellipsoid.

◆ geod_inverse()

Solve the inverse geodesic problem.

Parameters

g must have been initialized with a call to geod_init(). lat1 and lat2 should be in the range [−90°, 90°]. 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:

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);

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

Referenced by main().

◆ geod_geninverse()

The general inverse geodesic calculation.

Parameters

Returns

a12 arc length from point 1 to point 2 (degrees).

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

◆ geod_lineinit()

Initialize a geod_geodesicline object.

Parameters

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

The geod_mask values are:

• 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).

When initialized by this function, point 3 is undefined (l->s13 = l->a13 = NaN).

◆ geod_directline()

Initialize a geod_geodesicline object in terms of the direct geodesic problem.

Parameters

This function sets point 3 of the geod_geodesicline to correspond to point 2 of the direct geodesic problem. See geod_lineinit() for more information.

◆ geod_gendirectline()

Initialize a geod_geodesicline object in terms of the direct geodesic problem specified in terms of either distance or arc length.

Parameters

This function sets point 3 of the geod_geodesicline to correspond to point 2 of the direct geodesic problem. See geod_lineinit() for more information.

◆ geod_inverseline()

Initialize a geod_geodesicline object in terms of the inverse geodesic problem.

Parameters

This function sets point 3 of the geod_geodesicline to correspond to point 2 of the inverse geodesic problem. See geod_lineinit() for more information.

Referenced by main().

◆ geod_position()

Compute the position along a geod_geodesicline.

Parameters

l must have been initialized with a call, e.g., to geod_lineinit(), with caps |= GEOD_DISTANCE_IN (or caps = 0). 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():

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():

double lat[101], lon[101];

int i;

geod_init(&g, 6378137, 1/298.257223563);

for (i = 0; i <= 100; ++i) {

geod_position(&l, i * l.s13 * 0.01, lat + i, lon + i, 0);

printf("%.5f %.5f\n", lat[i], lon[i]);

}

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

void GEOD_DLL geod_inverseline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, unsigned caps)

◆ geod_genposition()

The general position function.

Parameters

Returns

a12 arc length from point 1 to point 2 (degrees).

l must have been initialized with a call to geod_lineinit() with caps |= GEOD_DISTANCE_IN. The value azi2 returned is 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.

With flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so that the quantity lon2 − lon1 indicates how many times and in what sense the geodesic encircles the ellipsoid.

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

double lat[101], lon[101];

int i;

geod_init(&g, 6378137, 1/298.257223563);

for (i = 0; i <= 100; ++i) {

lat + i, lon + i, 0, 0, 0, 0, 0, 0);

printf("%.5f %.5f\n", lat[i], lon[i]);

}

double GEOD_DLL geod_genposition(const struct geod_geodesicline *l, unsigned flags, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)

Referenced by main().

◆ geod_setdistance()

Specify position of point 3 in terms of distance.

Parameters

This is only useful if the geod_geodesicline object has been constructed with caps |= GEOD_DISTANCE_IN.

◆ geod_gensetdistance()

◆ geod_polygon_init()

Initialize a geod_polygon object.

Parameters

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().

The area and perimeter are accumulated at two times the standard floating point precision to guard against the loss of accuracy with many-sided polygons. At any point you can ask for the perimeter and area so far.

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

Referenced by main().

◆ geod_polygon_clear()

Clear the polygon, allowing a new polygon to be started.

Parameters

◆ geod_polygon_addpoint()

Add a point to the polygon or polyline.

Parameters

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°].

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

Referenced by main().

◆ geod_polygon_addedge()

Add an edge to the polygon or polyline.

Parameters

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. 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_compute()

Return the results for a polygon.

Parameters

Returns

the number of points.

The area and perimeter are accumulated at two times the standard floating point precision to guard against the loss of accuracy with many-sided polygons. Arbitrarily complex polygons are allowed. In the case of self-intersecting polygons the area is accumulated "algebraically", e.g., the areas of the 2 loops in a figure-8 polygon will partially cancel. 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.

More points can be added to the polygon after this call.

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;

geod_init(&g, 6378137, 1/298.257223563);

printf("%d %.8f %.3f\n", n, P, A);

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

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

void GEOD_DLL geod_polygon_init(struct geod_polygon *p, int polylinep)

Referenced by main().

◆ geod_polygon_testpoint()

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

Returns

the number of points.

lat should be in the range [−90°, 90°].

◆ geod_polygon_testedge()

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

Returns

the number of points.

◆ geod_polygonarea()

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

Parameters

lats should be in the range [−90°, 90°].

Arbitrarily complex polygons are allowed. In the case self-intersecting of polygons the area is accumulated "algebraically", e.g., the areas of the 2 loops in a figure-8 polygon will partially cancel. 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 Antarctica:

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};

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);

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