From 3175dce5fcec54fef6299cc8b595d84db5b29cf4 Mon Sep 17 00:00:00 2001 From: Aaron Puchert Date: Fri, 20 Oct 2017 12:58:10 +0200 Subject: Replace some asserts by compile-time assertions Decisions regarding the size of certain types are always made at compile time, hence this is where the assertion should take place. If the assertions fail, there is now a compiler error. --- src/PJ_cart.c | 1 - src/PJ_helmert.c | 1 - src/PJ_horner.c | 1 - src/PJ_pipeline.c | 1 - src/nad2bin.c | 12 ++++----- src/nad_init.c | 10 ------- src/pj_gridinfo.c | 78 +++++++++-------------------------------------------- src/pj_gridlist.c | 10 ------- src/proj_api.h | 2 +- src/proj_internal.h | 2 ++ src/projects.h | 11 +++++++- 11 files changed, 31 insertions(+), 98 deletions(-) diff --git a/src/PJ_cart.c b/src/PJ_cart.c index cd1995c1..12b5876a 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -43,7 +43,6 @@ #define PJ_LIB__ #include "proj_internal.h" #include "projects.h" -#include #include #include #include diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c index d0fd8222..6c4b520c 100644 --- a/src/PJ_helmert.c +++ b/src/PJ_helmert.c @@ -43,7 +43,6 @@ Last update: 2017-05-15 ***********************************************************************/ #define PJ_LIB__ -#include #include #include #include "proj_internal.h" diff --git a/src/PJ_horner.c b/src/PJ_horner.c index 6f1450cd..fe0452d6 100644 --- a/src/PJ_horner.c +++ b/src/PJ_horner.c @@ -78,7 +78,6 @@ #define PJ_LIB__ #include "proj_internal.h" #include "projects.h" -#include #include #include #include diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c index 35f79213..c8ce8582 100644 --- a/src/PJ_pipeline.c +++ b/src/PJ_pipeline.c @@ -102,7 +102,6 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-20 #include "proj_internal.h" #include "projects.h" -#include #include #include PROJ_HEAD(pipeline, "Transformation pipeline manager"); diff --git a/src/nad2bin.c b/src/nad2bin.c index 0001189c..8401df79 100644 --- a/src/nad2bin.c +++ b/src/nad2bin.c @@ -1,10 +1,10 @@ /* Convert bivariate ASCII NAD27 to NAD83 tables to NTv2 binary structure */ #include #include -#include #define PJ_LIB__ -#include +#include "proj_internal.h" +#include "projects.h" #define U_SEC_TO_RAD 4.848136811095359935899141023e-12 /************************************************************************/ @@ -177,8 +177,8 @@ int main(int argc, char **argv) { exit(2); } - assert( MAX_TAB_ID == 80 ); - assert( sizeof(int) == 4 ); /* for ct.lim.lam/phi */ + STATIC_ASSERT( MAX_TAB_ID == 80 ); + STATIC_ASSERT( sizeof(pj_int32) == 4 ); /* for ct.lim.lam/phi */ memset( header, 0, sizeof(header) ); @@ -267,13 +267,13 @@ int main(int argc, char **argv) { { unsigned char achHeader[11*16]; double dfValue; - int nGSCount = ct.lim.lam * ct.lim.phi; + pj_int32 nGSCount = ct.lim.lam * ct.lim.phi; LP ur; ur.lam = ct.ll.lam + (ct.lim.lam-1) * ct.del.lam; ur.phi = ct.ll.phi + (ct.lim.phi-1) * ct.del.phi; - assert( sizeof(nGSCount) == 4 ); + STATIC_ASSERT( sizeof(nGSCount) == 4 ); memset( achHeader, 0, sizeof(achHeader) ); diff --git a/src/nad_init.c b/src/nad_init.c index a9082f8f..99342aa5 100644 --- a/src/nad_init.c +++ b/src/nad_init.c @@ -32,16 +32,6 @@ #include #include -#ifdef _WIN32_WCE -/* assert.h includes all Windows API headers and causes 'LP' name clash. - * Here assert we disable assert() for Windows CE. - * TODO - mloskot: re-implement porting friendly assert - */ -# define assert(exp) ((void)0) -#else -# include -#endif /* _WIN32_WCE */ - /************************************************************************/ /* swap_words() */ /* */ diff --git a/src/pj_gridinfo.c b/src/pj_gridinfo.c index ad4695ca..9b9a8d82 100644 --- a/src/pj_gridinfo.c +++ b/src/pj_gridinfo.c @@ -28,21 +28,12 @@ #define PJ_LIB__ -#include +#include "proj_internal.h" +#include "projects.h" #include #include #include -#ifdef _WIN32_WCE -/* assert.h includes all Windows API headers and causes 'LP' name clash. - * Here assert we disable assert() for Windows CE. - * TODO - mloskot: re-implement porting friendly assert - */ -# define assert(exp) ((void)0) -#else -# include -#endif /* _WIN32_WCE */ - /************************************************************************/ /* swap_words() */ /* */ @@ -431,23 +422,8 @@ static int pj_gridinfo_init_ntv2( projCtx ctx, PAFile fid, PJ_GRIDINFO *gilist ) int num_subfiles, subfile; int must_swap; - assert( sizeof(int) == 4 ); - assert( sizeof(double) == 8 ); -#ifdef _MSC_VER -#pragma warning( push ) -/* disable conditional expression is constant */ -#pragma warning( disable : 4127 ) -#endif - if( sizeof(int) != 4 || sizeof(double) != 8 ) -#ifdef _MSC_VER -#pragma warning( pop ) -#endif - { - pj_log( ctx, PJ_LOG_ERROR, - "basic types of inappropraiate size in pj_gridinfo_init_ntv2()" ); - pj_ctx_set_errno( ctx, -38 ); - return 0; - } + STATIC_ASSERT( sizeof(pj_int32) == 4 ); + STATIC_ASSERT( sizeof(double) == 8 ); /* -------------------------------------------------------------------- */ /* Read the overview header. */ @@ -541,8 +517,8 @@ static int pj_gridinfo_init_ntv2( projCtx ctx, PAFile fid, PJ_GRIDINFO *gilist ) ct->del.lam = *((double *) (header+9*16+8)); ct->del.phi = *((double *) (header+8*16+8)); - ct->lim.lam = (int) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1; - ct->lim.phi = (int) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1; + ct->lim.lam = (pj_int32) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1; + ct->lim.phi = (pj_int32) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1; pj_log( ctx, PJ_LOG_DEBUG_MINOR, "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)\n", @@ -669,23 +645,8 @@ static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi ) struct CTABLE *ct; LP ur; - assert( sizeof(int) == 4 ); - assert( sizeof(double) == 8 ); -#ifdef _MSC_VER -#pragma warning( push ) -/* disable conditional expression is constant */ -#pragma warning( disable : 4127 ) -#endif - if( sizeof(int) != 4 || sizeof(double) != 8 ) -#ifdef _MSC_VER -#pragma warning( pop ) -#endif - { - pj_log( ctx, PJ_LOG_ERROR, - "basic types of inappropraiate size in nad_load_ntv1()" ); - pj_ctx_set_errno( ctx, -38 ); - return 0; - } + STATIC_ASSERT( sizeof(pj_int32) == 4 ); + STATIC_ASSERT( sizeof(double) == 8 ); /* -------------------------------------------------------------------- */ /* Read the header. */ @@ -734,8 +695,8 @@ static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi ) ur.phi = *((double *) (header+40)); ct->del.lam = *((double *) (header+104)); ct->del.phi = *((double *) (header+88)); - ct->lim.lam = (int) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1; - ct->lim.phi = (int) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1; + ct->lim.lam = (pj_int32) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1; + ct->lim.phi = (pj_int32) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1; pj_log( ctx, PJ_LOG_DEBUG_MINOR, "NTv1 %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", @@ -769,23 +730,8 @@ static int pj_gridinfo_init_gtx( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi ) double xorigin,yorigin,xstep,ystep; int rows, columns; - assert( sizeof(int) == 4 ); - assert( sizeof(double) == 8 ); -#ifdef _MSC_VER -#pragma warning( push ) -/* disable conditional expression is constant */ -#pragma warning( disable : 4127 ) -#endif - if( sizeof(int) != 4 || sizeof(double) != 8 ) -#ifdef _MSC_VER -#pragma warning( pop ) -#endif - { - pj_log( ctx, PJ_LOG_ERROR, - "basic types of inappropraiate size in nad_load_gtx()" ); - pj_ctx_set_errno( ctx, -38 ); - return 0; - } + STATIC_ASSERT( sizeof(pj_int32) == 4 ); + STATIC_ASSERT( sizeof(double) == 8 ); /* -------------------------------------------------------------------- */ /* Read the header. */ diff --git a/src/pj_gridlist.c b/src/pj_gridlist.c index 4617591c..fbefdcea 100644 --- a/src/pj_gridlist.c +++ b/src/pj_gridlist.c @@ -33,16 +33,6 @@ #include #include -#ifdef _WIN32_WCE -/* assert.h includes all Windows API headers and causes 'LP' name clash. - * Here assert we disable assert() for Windows CE. - * TODO - mloskot: re-implement porting friendly assert - */ -# define assert(exp) ((void)0) -#else -# include -#endif /* _WIN32_WCE */ - static PJ_GRIDINFO *grid_list = NULL; #define PJ_MAX_PATH_LENGTH 1024 diff --git a/src/proj_api.h b/src/proj_api.h index b57d9f38..e5209615 100644 --- a/src/proj_api.h +++ b/src/proj_api.h @@ -69,10 +69,10 @@ extern "C" { -#ifndef PROJ_H extern char const pj_release[]; /* global release id string */ extern int pj_errno; /* global error return code */ +#ifndef PROJ_H /* In proj.h these macros are replaced by the enumeration pj_log_level */ #define PJ_LOG_NONE 0 #define PJ_LOG_ERROR 1 diff --git a/src/proj_internal.h b/src/proj_internal.h index fd6dc75d..4e70e690 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -48,6 +48,8 @@ extern "C" { #endif +#define STATIC_ASSERT(COND) ((void)sizeof(char[(COND) ? 1 : -1])) + #ifndef PJ_TODEG #define PJ_TODEG(rad) ((rad)*180.0/M_PI) #endif diff --git a/src/projects.h b/src/projects.h index 99df0e3e..8c0f81fa 100644 --- a/src/projects.h +++ b/src/projects.h @@ -43,6 +43,7 @@ #endif /* standard inclusions */ +#include #include #include #include @@ -78,6 +79,14 @@ extern "C" { # define ABS(x) ((x<0) ? (-1*(x)) : x) #endif +#if INT_MAX == 2147483647 +typedef int pj_int32; +#elif LONG_MAX == 2147483647 +typedef long pj_int32; +#else +#warning It seems no 32-bit integer type is available +#endif + /* maximum path/filename */ #ifndef MAX_PATH_FILENAME #define MAX_PATH_FILENAME 1024 @@ -618,7 +627,7 @@ PJ *pj_projection_specific_setup_##name (PJ *P) #define MAX_TAB_ID 80 typedef struct { float lam, phi; } FLP; -typedef struct { int lam, phi; } ILP; +typedef struct { pj_int32 lam, phi; } ILP; struct CTABLE { char id[MAX_TAB_ID]; /* ascii info */ -- cgit v1.2.3 From a141977eb6a476b4dfc7e35aaa26531f01f2fff9 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Mon, 27 Nov 2017 11:36:31 -0500 Subject: Update documentation on geodesics replacing an ancient page on geodesic calculations with something more helpful to the beginning user. --- docs/source/apps/geod.rst | 10 +- docs/source/faq.rst | 2 +- docs/source/geodesic.rst | 424 +++++++++++++++++++++++----------------------- man/man3/geodesic.3 | 6 +- 4 files changed, 221 insertions(+), 221 deletions(-) diff --git a/docs/source/apps/geod.rst b/docs/source/apps/geod.rst index 4620a6cc..49fa4925 100644 --- a/docs/source/apps/geod.rst +++ b/docs/source/apps/geod.rst @@ -163,11 +163,9 @@ which gives: Further reading *************** -#. `GeographicLib `_ - -#. `C. F. F. Karney, Algorithms for Geodesics, J. Geodesy 87, 43-55 (2013) `_. - `Addendum `_ - -#. `The online geodesic bibliography `_ +#. `GeographicLib `_. +#. C. F. F. Karney, `Algorithms for Geodesics `_, J. Geodesy **87**\ (1), 43–55 (2013); + `addenda `_. +#. `A geodesic bibliography `_. diff --git a/docs/source/faq.rst b/docs/source/faq.rst index 0240df1d..14ceda85 100644 --- a/docs/source/faq.rst +++ b/docs/source/faq.rst @@ -237,7 +237,7 @@ How do I calculate distances/directions on the surface of the earth? -------------------------------------------------------------------------------- These are called geodesic calculations. There is a page about it here: -[wiki:GeodesicCalculations] +:ref:`geodesic`. What is the HEALPix projection and how can I use it? -------------------------------------------------------------------------------- diff --git a/docs/source/geodesic.rst b/docs/source/geodesic.rst index fac3649c..29327ae8 100644 --- a/docs/source/geodesic.rst +++ b/docs/source/geodesic.rst @@ -1,219 +1,217 @@ .. _geodesic: -================================================================================ -Geodesic Calculations -================================================================================ +Geodesic calculations +===================== .. contents:: Contents - :depth: 3 + :depth: 2 :backlinks: none - -Geodesic Calculations --------------------------------------------------------------------------------- - -Geodesic calculations are calculations along lines (great circle) on the -surface of the earth. They can answer questions like: - - * What is the distance between these two points? - * If I travel X meters from point A at bearing phi, where will I be. They are - done in native lat-long coordinates, rather than in projected coordinates. - -Relevant mailing list threads -................................................................................ - - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3361 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3375 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3435 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3588 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3925 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/4047 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/4083 - -Terminology --------------------------------------------------------------------------------- - -The shortest distance on the surface of a solid is generally termed a geodesic, -be it an ellipsoid of revolution, aposphere, etc. On a sphere, the geodesic is -termed a Great Circle. - -HOWEVER, when computing the distance between two points using a projected -coordinate system, that is a conformal projection such as Transverse Mercator, -Oblique Mercator, Normal Mercator, Stereographic, or Lambert Conformal Conic - -that then is a GRID distance which can be converted to an equivalent GEODETIC -distance using the function for "Scale Factor at a Point." The conversion is -then termed "Grid Distance to Geodetic Distance," even though it will not be as -exactly correct as a true ellipsoidal geodesic. Closer to the truth with a TM -than with a Lambert or other conformal projection, but still not exactly "on." - - -So, it can be termed "geodetic distance" or a "geodesic distance," depending -on just how you got there ... - - -The Math --------------------------------------------------------------------------------- - -Spherical Approximation -................................................................................ - -The simplest way to compute geodesics is using a sphere as an approximation for -the earth. This from Mikael Rittri on the Proj mailing list: - - If 1 percent accuracy is enough, I think you can use spherical formulas - with a fixed Earth radius. You can find good formulas in the Aviation - Formulary of Ed Williams, http://williams.best.vwh.net/avform.htm. - - For the fixed Earth radius, I would choose the average of the: - - :math:`c` = radius of curvature at the poles, - - :math:`\frac{b^2}{a}` = radius of curvature in a meridian plane at the equator, - - since these are the extreme values for the local radius of curvature of the - earth ellipsoid. - - If your coordinates are given in WGS84, then - - :math:`c` = 6 399 593.626 m, - - :math:`\frac{b^2}{a}` = 6 335 439.327 m, - - (see http://home.online.no/~sigurdhu/WGS84_Eng.html) so their average is 6,367,516.477 m. - The maximal error for distance calculation should then be less than 0.51 percent. - - When computing the azimuth between two points by the spherical formulas, I - think the maximal error on WGS84 will be 0.2 degrees, at least if the - points are not too far away (less than 1000 km apart, say). The error - should be maximal near the equator, for azimuths near northeast etc. - -I am not sure about the spherical errors for the forward geodetic problem: -point positioning given initial point, distance and azimuth. - -Ellipsoidal Approximation -................................................................................ - -For more accuracy, the earth can be approximated with an ellipsoid, -complicating the math somewhat. See the wikipedia page, `Geodesics on an -ellipsoid `__, for -more information. - -Thaddeus Vincenty's method, April 1975 -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -For a very good procedure to calculate inter point distances see: - -http://www.ngs.noaa.gov/PC_PROD/Inv_Fwd/ (Fortran code, DOS executables, and an online app) - -and algorithm details published in: `Vincenty, T. (1975) `__ - -Javascript code -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Chris Veness has coded Vincenty's formulas as !JavaScript. - -distance: http://www.movable-type.co.uk/scripts/latlong-vincenty.html - -direct: http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html - -C code -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -From Gerald Evenden: a library of the converted NGS Vincenty geodesic procedure -and an application program, 'geodesic'. In the case of a spherical earth -[Snyder1987]_'s preferred equations are used. - -* http://article.gmane.org/gmane.comp.gis.proj-4.devel/3588/ - -The link in this message is broken. The correct URL is -http://home.comcast.net/~gevenden56/proj/ - -Earlier Mr. Evenden had posted to the PROJ.4 mailing list this code for -determination of true distance and respective forward and back azimuths between -two points on the ellipsoid. Good for any pair of points that are not -antipodal. -Later he posted that this was not in fact the translation of NGS FORTRAN code, -but something else. But, for what it's worth, here is the posted code (source -unknown): - -* http://article.gmane.org/gmane.comp.gis.proj-4.devel/3478 - - -PROJ.4 - geod program -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - -The PROJ.4 [wiki:man_geod geod] program can be used for great circle distances -on an ellipsoid. As of proj version 4.9.0, this uses a translation of -GeographicLib::Geodesic (see below) into C. The underlying geodesic -calculation API is exposed as part of the PROJ.4 library (via the geodesic.h -header). Prior to version 4.9.0, the algorithm documented here was used: -` -Paul D. Thomas, 1970 -Spheroidal Geodesics, Reference Systems, and Local Geometry" -U.S. Naval Oceanographic Office, p. 162 -Engineering Library 526.3 T36s - -http://handle.dtic.mil/100.2/AD0703541 - -GeographicLib::Geodesic -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -Charles Karney has written a C++ class to do geodesic calculations and a -utility GeodSolve to call it. See - -* http://geographiclib.sourceforge.net/geod.html - -An online version of GeodSolve is available at - -* http://geographiclib.sourceforge.net/cgi-bin/GeodSolve - -This is an attempt to do geodesic calculations "right", i.e., - -* accurate to round-off (i.e., about 15 nm); -* inverse solution always succeeds (even for near anti-podal points); -* reasonably fast (comparable in speed to Vincenty); -* differential properties of geodesics are computed (these give the scales of - geodesic projections); -* the area between a geodesic and the equator is computed (allowing the - area of geodesic polygons to be found); -* included also is an implementation in terms of elliptic integrals which - can deal with ellipsoids with 0.01 < b/a < 100. - -A JavaScript implementation is included, see - -* `geo-calc `__, - a text interface to geodesic calculations; -* `geod-google `__, - a tool for drawing geodesics on Google Maps. - -Implementations in `Python `__, -`Matlab `__, -`C `__, -`Fortran `__ , and -`Java `__ are also available. - -The algorithms are described in - * C. F. F. Karney, `Algorithms for gedesics `__, - J. Geodesy '''87'''(1), 43-55 (2013), - DOI: `10.1007/s00190-012-0578-z `__; `geo-addenda.html `__. - -Triaxial Ellipsoid -................................................................................ - -A triaxial ellipsoid is a marginally better approximation to the shape of the earth -than an ellipsoid of revolution. -The problem of geodesics on a triaxial ellipsoid was solved by Jacobi in 1838. -For a discussion of this problem see -* http://geographiclib.sourceforge.net/html/triaxial.html -* the wikipedia entry: `Geodesics on a triaxial ellipsoid `__ - -The History --------------------------------------------------------------------------------- - -The bibliography of papers on the geodesic problem for an ellipsoid is -available at - -* http://geographiclib.sourceforge.net/geodesic-papers/biblio.html - -this includes links to online copies of the papers. +Introduction +------------ + +Consider a ellipsoid of revolution with equatorial radius *a*, polar +semi-axis *b*, and flattening *f* = (*a* − *b*)/*a* . Points on +the surface of the ellipsoid are characterized by their latitude φ +and longitude λ. (Note that latitude here means the +*geographical latitude*, the angle between the normal to the ellipsoid +and the equatorial plane). + +The shortest path between two points on the ellipsoid at +(φ\ :sub:`1`, λ\ :sub:`1`) and (φ\ :sub:`2`, λ\ :sub:`2`) +is called the geodesic. Its length is +*s*\ :sub:`12` and the geodesic from point 1 to point 2 has forward +azimuths α\ :sub:`1` and α\ :sub:`2` at the two end +points. In this figure, we have λ\ :sub:`12` = +λ\ :sub:`2` − λ\ :sub:`1`. + + .. raw:: html + +
+ Figure from wikipedia +
+ +A geodesic can be extended indefinitely by requiring that any +sufficiently small segment is a shortest path; geodesics are also the +straightest curves on the surface. + +Solution of geodesic programs +----------------------------- + +Traditionally two geodesic problems are considered: + +* the direct problem — given φ\ :sub:`1`, + λ\ :sub:`1`, α\ :sub:`1`, *s*\ :sub:`12`, + determine φ\ :sub:`2`, λ\ :sub:`2`, and + α\ :sub:`2`. + +* the inverse problem — given φ\ :sub:`1`, + λ\ :sub:`1`, φ\ :sub:`2`, λ\ :sub:`2`, + determine *s*\ :sub:`12`, α\ :sub:`1`, and + α\ :sub:`2`. + +PROJ incorporates `C library for Geodesics +`_ from `GeographicLib +`_. This library provides +routines to solve the direct and inverse geodesic problems. Full double +precision accuracy is maintained provided that −0.02 < *f* < 0.02. Refer +to the + + `application programming interface + `_ + +for full documentation. A brief summary of the routines is given in +geodesic(3). + +The interface to the geodesic routines differ in two respects from the +rest of PROJ: + +* angles (latitudes, longitudes, and azimuths) are in degrees (instead + of in radians); +* the shape of ellipsoid is specified by the flattening *f*; this can + be negative to denote a prolate ellipsoid; setting *f* = 0 corresponds + to a sphere in which case the geodesic becomes a great circle. + +PROJ also includes a command line tool, :ref:`geod`\ (1), for performing +simple geodesic calculations. + +Additional properties +--------------------- + +The routines also calculate several other quantities of interest + +* *S*\ :sub:`12` 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 + (φ\ :sub:`1`,λ\ :sub:`1`), (0,λ\ :sub:`1`), + (0,λ\ :sub:`2`), and + (φ\ :sub:`2`,λ\ :sub:`2`). It is given in + meters\ :sup:`2`. +* *m*\ :sub:`12`, the reduced length of the geodesic is defined such + that if the initial azimuth is perturbed by *d*\ α\ :sub:`1` + (radians) then the second point is displaced by *m*\ :sub:`12` + *d*\ α\ :sub:`1` in the direction perpendicular to the + geodesic. *m*\ :sub:`12` is given in meters. On a curved surface + the reduced length obeys a symmetry relation, *m*\ :sub:`12` + + *m*\ :sub:`21` = 0. On a flat surface, we have *m*\ :sub:`12` = + *s*\ :sub:`12`. +* *M*\ :sub:`12` and *M*\ :sub:`21` 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 *M*\ :sub:`12` *dt* at + point 2. *M*\ :sub:`21` is defined similarly (with the geodesics + being parallel to one another at point 2). *M*\ :sub:`12` and + *M*\ :sub:`21` are dimensionless quantities. On a flat surface, + we have *M*\ :sub:`12` = *M*\ :sub:`21` = 1. +* σ\ :sub:`12` is the arc length on the auxiliary sphere. + This is a construct for converting the problem to one in spherical + trigonometry. 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: + +* *s*\ :sub:`13` = *s*\ :sub:`12` + *s*\ :sub:`23` +* σ\ :sub:`13` = σ\ :sub:`12` + σ\ :sub:`23` +* *S*\ :sub:`13` = *S*\ :sub:`12` + *S*\ :sub:`23` +* *m*\ :sub:`13` = *m*\ :sub:`12`\ *M*\ :sub:`23` + + *m*\ :sub:`23`\ *M*\ :sub:`21` +* *M*\ :sub:`13` = *M*\ :sub:`12`\ *M*\ :sub:`23` − + (1 − *M*\ :sub:`12`\ *M*\ :sub:`21`) + *m*\ :sub:`23`/*m*\ :sub:`12` +* *M*\ :sub:`31` = *M*\ :sub:`32`\ *M*\ :sub:`21` − + (1 − *M*\ :sub:`23`\ *M*\ :sub:`32`) + *m*\ :sub:`12`/*m*\ :sub:`23` + +Multiple shortest geodesics +--------------------------- + +The shortest distance found by solving 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: + +* φ\ :sub:`1` = −φ\ :sub:`2` (with neither point at + a pole). If α\ :sub:`1` = α\ :sub:`2`, the geodesic + is unique. Otherwise there are two geodesics and the second one is + obtained by setting + [α\ :sub:`1`,α\ :sub:`2`] ← [α\ :sub:`2`,α\ :sub:`1`], + [*M*\ :sub:`12`,\ *M*\ :sub:`21`] ← [*M*\ :sub:`21`,\ *M*\ :sub:`12`], + *S*\ :sub:`12` ← −\ *S*\ :sub:`12`. + (This occurs when the longitude difference is near ±180° for oblate + ellipsoids.) +* λ\ :sub:`2` = λ\ :sub:`1` ± 180° (with + neither point at a pole). If α\ :sub:`1` = 0° or + ±180°, the geodesic is unique. Otherwise there are two + geodesics and the second one is obtained by setting + [α\ :sub:`1`,α\ :sub:`2`] ← [−α\ :sub:`1`,−α\ :sub:`2`], + *S*\ :sub:`12` ← −\ *S*\ :sub:`12`. (This occurs when + φ\ :sub:`2` is near −φ\ :sub:`1` for prolate + ellipsoids.) +* Points 1 and 2 at opposite poles. There are infinitely many + geodesics which can be generated by setting + [α\ :sub:`1`,α\ :sub:`2`] ← + [α\ :sub:`1`,α\ :sub:`2`] + [δ,−δ], for arbitrary δ. + (For spheres, this prescription applies when points 1 and 2 are + antipodal.) +* *s*\ :sub:`12` = 0 (coincident points). There are infinitely many + geodesics which can be generated by setting + [α\ :sub:`1`,α\ :sub:`2`]_← + [α\ :sub:`1`,α\ :sub:`2`]_+ [δ,δ], for + arbitrary δ. + +Background +---------- + +The algorithms implemented by this package are given in Karney (2013) +and are based on Bessel (1825) and Helmert (1880); the algorithm for +areas is based on Danielsen (1989). These improve on the work of +Vincenty (1975) in the following respects: + +* The results are accurate to round-off for terrestrial ellipsoids (the + error in the distance is less then 15 nanometers, compared to 0.1 mm + for Vincenty). +* The solution of the inverse problem is always found. (Vincenty's + method fails to converge for nearly antipodal points.) +* The routines calculate differential and integral properties of a + geodesic. This allows, for example, the area of a geodesic polygon to + be computed. + +References +---------- + +* F. W. Bessel, + `The calculation of longitude and latitude from geodesic measurements (1825) + `_, + Astron. Nachr. **331**\ (8), 852–861 (2010), + translated by C. F. F. Karney and R. E. Deakin. +* F. R. Helmert, + `Mathematical and Physical Theories of Higher Geodesy, Vol 1 + `_, + (Teubner, Leipzig, 1880), Chaps. 5–7. +* T. Vincenty, + `Direct and inverse solutions of geodesics on the ellipsoid with + application of nested equations + `_, + Survey Review **23**\ (176), 88–93 (1975). +* J. Danielsen, + `The area under the geodesic + `_, + Survey Review **30**\ (232), 61–66 (1989). +* C. F. F. Karney, + `Algorithms for geodesics + `_, + J. Geodesy **87**\ (1) 43–55 (2013); + `addenda `_. +* C. F. F. Karney, + `Geodesics on an ellipsoid of revolution + `_, + Feb. 2011; + `errata + `_. +* `A geodesic bibliography + `_. +* The wikipedia page, + `Geodesics on an ellipsoid + `_. diff --git a/man/man3/geodesic.3 b/man/man3/geodesic.3 index 938eed68..52f9b042 100644 --- a/man/man3/geodesic.3 +++ b/man/man3/geodesic.3 @@ -53,7 +53,9 @@ measure angles (latitudes, longitudes, and azimuths) in degrees, unlike the rest of the \fBproj\fR library, which uses radians. The documentation for this library is included in geodesic.h. A formatted version of the documentation is available at -https://geographiclib.sourceforge.io/1.49/C +https://geographiclib.sourceforge.io/1.49/C. Detailed documentation of +the interface is given at +https://geographiclib.sourceforge.io/1.49/C/geodesic_8h.html. .SH EXAMPLE The following program reads in lines with the coordinates for two points in decimal degrees (\fIlat1\fR, \fIlon1\fR, \fIlat2\fR, \fIlon2\fR) and @@ -88,6 +90,8 @@ libproj.a \- library of projections and support procedures Full online documentation for \fBgeodesic(3)\fR, .br https://geographiclib.sourceforge.io/1.49/C +.br +https://geographiclib.sourceforge.io/1.49/C/geodesic_8h.html .PP .B geod(1) .PP -- cgit v1.2.3 From 8c6fe8b673c9876e7301d4a890403a931540c17c Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Tue, 28 Nov 2017 16:09:17 +1300 Subject: Add tests for GDA2020 (#688) Add GDA2020 tests, and correct the PJ_helmert bug they revealed --- appveyor.yml | 1 + src/PJ_helmert.c | 3 ++- test/gie/GDA.gie | 69 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ travis/install.sh | 1 + 4 files changed, 73 insertions(+), 1 deletion(-) create mode 100644 test/gie/GDA.gie diff --git a/appveyor.yml b/appveyor.yml index 43669a40..f8bbd7ef 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -63,6 +63,7 @@ test_script: - gie.exe ..\test\gie\deformation.gie - gie.exe ..\test\gie\axisswap.gie - gie.exe ..\test\gie\ellipsoid.gie + - gie.exe ..\test\gie\GDA.gie deploy: off diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c index d0fd8222..79c7ad53 100644 --- a/src/PJ_helmert.c +++ b/src/PJ_helmert.c @@ -572,7 +572,8 @@ PJ *TRANSFORMATION(helmert, 0) { proj_log_debug(P, "ds=% 3.5f epoch=% 5.5f tobs=% 5.5f", Q->dscale, Q->epoch, Q->t_obs); } - if ((Q->opk.o==0) && (Q->opk.p==0) && (Q->opk.k==0) && (Q->scale==0)) { + if ((Q->opk.o==0) && (Q->opk.p==0) && (Q->opk.k==0) && (Q->scale==0) && + (Q->dopk.o==0) && (Q->dopk.p==0) && (Q->dopk.k==0)) { Q->no_rotation = 1; return P; } diff --git a/test/gie/GDA.gie b/test/gie/GDA.gie new file mode 100644 index 00000000..7ccf6b58 --- /dev/null +++ b/test/gie/GDA.gie @@ -0,0 +1,69 @@ +----------------------------------------------------------------------------------- +Australian datum transformations +----------------------------------------------------------------------------------- +Based on material from: + +Intergovernmental Committee on Surveying and Mapping (ICSM) +Permanent Committee on Geodesy (PCG): + +Geocentric Datum of Australia 2020 Technical Manual +Version 1.0, 25 July 2017 + +Which is distributed under Creative Commons CC-BY 4.0 + +These tests will probably be useful as a template for an AU setup file, defining +transformations for Australian systems, but I'm reluctant to provide such a file +directly - I believe it should come from official AU sources. + +Thomas Knudsen, thokn@sdfe.dk, 2017-11-27 +----------------------------------------------------------------------------------- + +BEGIN +----------------------------------------------------------------------------------- +GDA94 to GDA2020 +----------------------------------------------------------------------------------- +Just the Helmert transformation, to verify that we are within 100 um +----------------------------------------------------------------------------------- +operation proj=helmert x=0.06155 y=-0.01087 z=-0.04019 s=-0.009994 \ + rx=-0.0394924 ry=-0.0327221 rz=-0.0328979 +----------------------------------------------------------------------------------- +tolerance 75 um +accept -4052051.7643 4212836.2017 -2545106.0245 +expect -4052052.7379 4212835.9897 -2545104.5898 + + +----------------------------------------------------------------------------------- +GDA94 to GDA2020 +----------------------------------------------------------------------------------- +All the way from geographic-to-geographic +----------------------------------------------------------------------------------- +operation proj=pipeline ellps=GRS80 \ + step proj=cart \ + step proj=helmert x=0.06155 y=-0.01087 z=-0.04019 s=-0.009994 \ + rx=-0.0394924 ry=-0.0327221 rz=-0.0328979 \ + step proj=cart inv +----------------------------------------------------------------------------------- +accept 133.88551329 -23.67012389 603.3466 0 # Alice Springs GDA94 +expect 133.8855216 -23.67011014 603.2489 0 # Alice Springs GDA2020 + + +----------------------------------------------------------------------------------- +ITRF2014@2018 to GDA2020 +----------------------------------------------------------------------------------- +Just the Helmert transformation, to verify that we are within 100 um +----------------------------------------------------------------------------------- +operation proj=helmert x=0 y=0 z=0 s=0 rx=0 ry=0 rz=0 \ + dx=0 dy=0 dz=0 \ + drx=0.00150379 dry=0.00118346 drz=0.00120716 \ + ds=0 \ + epoch=2020.0 +----------------------------------------------------------------------------------- +tolerance 50 um + + # Alice Springs ITRF2014@2018.0 +accept -4052052.6588 4212835.9938 -2545104.6946 2018.0 + + # Alice Springs GDA2020 +expect -4052052.7373 4212835.9835 -2545104.5867 + +END diff --git a/travis/install.sh b/travis/install.sh index f71108c6..78da24de 100755 --- a/travis/install.sh +++ b/travis/install.sh @@ -85,6 +85,7 @@ PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/more_builtins.gie PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/deformation.gie PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/axisswap.gie PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/ellipsoid.gie +PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/GDA.gie # install & run the working GIGS test # create locations that pyproj understands -- cgit v1.2.3 From 1c1d04b45d76366f54e104f9346879fd48bfde8e Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Tue, 28 Nov 2017 18:34:02 +1300 Subject: Enable selective debug printout through negative PROJ_DEBUG values (#689) * Enable selective debug printout through negative PROJ_DEBUG values --- src/pj_ctx.c | 2 +- src/pj_internal.c | 2 +- src/pj_log.c | 11 ++++++++++- src/proj_api.h | 5 +++-- 4 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/pj_ctx.c b/src/pj_ctx.c index 326249ac..fc52f300 100644 --- a/src/pj_ctx.c +++ b/src/pj_ctx.c @@ -84,7 +84,7 @@ projCtx pj_get_default_ctx() if( getenv("PROJ_DEBUG") != NULL ) { - if( atoi(getenv("PROJ_DEBUG")) > 0 ) + if( atoi(getenv("PROJ_DEBUG")) >= -PJ_LOG_DEBUG_MINOR ) default_context.debug_level = atoi(getenv("PROJ_DEBUG")); else default_context.debug_level = PJ_LOG_DEBUG_MINOR; diff --git a/src/pj_internal.c b/src/pj_internal.c index efa9bd1f..8a5d2d15 100644 --- a/src/pj_internal.c +++ b/src/pj_internal.c @@ -183,7 +183,7 @@ enum proj_log_level proj_log_level (PJ_CONTEXT *ctx, enum proj_log_level log_lev ctx = pj_get_default_ctx(); if (0==ctx) return PJ_LOG_TELL; - previous = ctx->debug_level; + previous = abs (ctx->debug_level); if (PJ_LOG_TELL==log_level) return previous; ctx->debug_level = log_level; diff --git a/src/pj_log.c b/src/pj_log.c index 2525d050..bab5f078 100644 --- a/src/pj_log.c +++ b/src/pj_log.c @@ -51,8 +51,17 @@ void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args ) { char *msg_buf; + int debug_level = ctx->debug_level; + int shutup_unless_errno_set = debug_level < 0; - if( level > ctx->debug_level ) + /* For negative debug levels, we first start logging when errno is set */ + if (ctx->last_errno==0 && shutup_unless_errno_set) + return; + + if (debug_level < 0) + debug_level = -debug_level; + + if( level > debug_level ) return; msg_buf = (char *) malloc(100000); diff --git a/src/proj_api.h b/src/proj_api.h index b57d9f38..50999b71 100644 --- a/src/proj_api.h +++ b/src/proj_api.h @@ -72,15 +72,16 @@ extern "C" { #ifndef PROJ_H extern char const pj_release[]; /* global release id string */ extern int pj_errno; /* global error return code */ +#endif -/* In proj.h these macros are replaced by the enumeration pj_log_level */ +#ifndef PROJ_INTERNAL_H +/* replaced by enum proj_log_level in proj_internal.h */ #define PJ_LOG_NONE 0 #define PJ_LOG_ERROR 1 #define PJ_LOG_DEBUG_MAJOR 2 #define PJ_LOG_DEBUG_MINOR 3 #endif - #ifdef PROJ_API_H_NOT_INVOKED_AS_PRIMARY_API /* These make the function declarations below conform with classic proj */ typedef PJ *projPJ; /* projPJ is a pointer to PJ */ -- cgit v1.2.3 From b1ec195f96e1b88de59e2b8d7e36b8e52730392a Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Tue, 28 Nov 2017 07:06:29 +0100 Subject: use full error code checking on failure in more_builtins.gie --- test/gie/more_builtins.gie | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie index 8ac8a667..ea37a763 100644 --- a/test/gie/more_builtins.gie +++ b/test/gie/more_builtins.gie @@ -140,11 +140,11 @@ roundtrip 100 1 nm ------------------------------------------------------------------------------- Fail on purpose: +grids parameter is mandatory operation proj=vgridshift -expect failure +expect failure errno no_args Fail on purpose: open non-existing grid operation proj=vgridshift grids=nonexistinggrid.gtx -expect failure +expect failure errno failed_to_load_grid @@ -166,11 +166,11 @@ expect 173 -45 0 0 ------------------------------------------------------------------------------- Fail on purpose: +grids parameter is mandatory: operation proj=hgridshift -expect failure +expect failure errno no_args Fail on purpose: open non-existing grid: operation proj=hgridshift grids=@nonexistinggrid.gsb,anothernonexistinggrid.gsb -expect failure +expect failure errno failed_to_load_grid ------------------------------------------------------------------------------- @@ -277,7 +277,7 @@ expect 12 89.999999999989996 0 0 some less used options ------------------------------------------------------------------------------- operation proj=utm ellps=GRS80 zone=32 to_meter=0 -expect failure +expect failure errno unit_factor_less_than_0 operation proj=utm ellps=GRS80 zone=32 to_meter=10 accept 12 55 -- cgit v1.2.3 From e51a6c255822d0a1b1320a2f519610f8802a5e27 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Tue, 28 Nov 2017 08:16:42 +0100 Subject: Free formatting of PROJ key=value pairs --- src/proj_4D_api.c | 44 ++++++++++++++++++++++++++++++-------------- test/gie/GDA.gie | 15 ++++++++------- 2 files changed, 38 insertions(+), 21 deletions(-) diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c index 9a576997..5f4bf334 100644 --- a/src/proj_4D_api.c +++ b/src/proj_4D_api.c @@ -381,40 +381,60 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) { Create a new PJ object in the context ctx, using the given definition. If ctx==0, the default context is used, if definition==0, or invalid, a null-pointer is returned. The definition may use '+' as argument start indicator, as in - "+proj=utm +zone=32", or leave it out, as in "proj=utm zone=32" + "+proj=utm +zone=32", or leave it out, as in "proj=utm zone=32". + + It may even use free formatting "proj = utm; zone =32 ellps= GRS80". + Note that the semicolon separator is allowed, but not required. **************************************************************************************/ PJ *P; char *args, **argv; - int argc, i, j, n; + int argc, i, j, last, n; if (0==ctx) ctx = pj_get_default_ctx (); - /* make a copy that we can manipulate */ + /* Make a copy that we can manipulate */ n = (int) strlen (definition); args = (char *) malloc (n + 1); if (0==args) return 0; strcpy (args, definition); - /* all-in-one: count args, eliminate superfluous whitespace, 0-terminate substrings */ - for (i = j = argc = 0; i < n; ) { - /* skip prefix whitespace */ + /* All-in-one: count args, eliminate superfluous whitespace, 0-terminate substrings */ + for (i = j = argc = last = 0; i < n; ) { + + /* Skip prefix whitespace */ while (isspace (args[i])) i++; - /* skip at most one prefix '+' */ + /* Skip at most one prefix '+' */ if ('+'==args[i]) i++; - /* whitespace after a '+' is a syntax error - but by Postel's prescription, we ignore and go on */ + /* Whitespace after a '+' is a syntax error - but by Postel's prescription, we ignore and go on */ if (isspace (args[i])) continue; - /* move a whitespace delimited text string to the left, skipping over superfluous whitespace */ - while ((0!=args[i]) && (!isspace (args[i]))) + /* Move a whitespace delimited text string to the left, skipping over superfluous whitespace */ + while ((0!=args[i]) && (!isspace (args[i])) && (';'!=args[i])) args[j++] = args[i++]; + /* Skip postfix whitespace */ + while (isspace (args[i]) || ';'==args[i]) + i++; + + /* Greedy assignment operator: turn "a = b" into "a=b" */ + if ('='==args[i]) { + args[j++] = '='; + i++; + while (isspace (args[i])) + i++; + while ((0!=args[i]) && (!isspace (args[i])) && (';'!=args[i])) + args[j++] = args[i++]; + while (isspace (args[i]) || ';'==args[i]) + i++; + } + /* terminate string - if that makes j pass i (often the case for first arg), let i catch up */ args[j++] = 0; if (i < j) @@ -422,10 +442,6 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) { /* we finished another arg */ argc++; - - /* skip postfix whitespace */ - while (isspace (args[i])) - i++; } /* turn the massaged input into an array of strings */ diff --git a/test/gie/GDA.gie b/test/gie/GDA.gie index 7ccf6b58..3fe0f3e3 100644 --- a/test/gie/GDA.gie +++ b/test/gie/GDA.gie @@ -18,7 +18,7 @@ directly - I believe it should come from official AU sources. Thomas Knudsen, thokn@sdfe.dk, 2017-11-27 ----------------------------------------------------------------------------------- -BEGIN + ----------------------------------------------------------------------------------- GDA94 to GDA2020 ----------------------------------------------------------------------------------- @@ -52,11 +52,12 @@ ITRF2014@2018 to GDA2020 ----------------------------------------------------------------------------------- Just the Helmert transformation, to verify that we are within 100 um ----------------------------------------------------------------------------------- -operation proj=helmert x=0 y=0 z=0 s=0 rx=0 ry=0 rz=0 \ - dx=0 dy=0 dz=0 \ - drx=0.00150379 dry=0.00118346 drz=0.00120716 \ - ds=0 \ - epoch=2020.0 +operation proj = helmert; \ + x = 0; rx = 0; dx = 0; drx = 0.00150379; \ + y = 0; ry = 0; dy = 0; dry = 0.00118346; \ + z = 0; rz = 0; dz = 0; drz = 0.00120716; \ + \ + ds = 0; epoch=2020.0; ----------------------------------------------------------------------------------- tolerance 50 um @@ -66,4 +67,4 @@ accept -4052052.6588 4212835.9938 -2545104.6946 2018.0 # Alice Springs GDA2020 expect -4052052.7373 4212835.9835 -2545104.5867 -END + -- cgit v1.2.3 From 38f717fb6b9838fd30ef1388bc66cd14005309bd Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Wed, 29 Nov 2017 20:25:40 +0100 Subject: Fix bad typedef that surfaces when including proj.h and proj_api.h in the same file (#687) Fix bad typedef that surfaces when including proj.h and proj_api.h in the same file --- src/proj_api.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/proj_api.h b/src/proj_api.h index f5ced282..597a2589 100644 --- a/src/proj_api.h +++ b/src/proj_api.h @@ -83,7 +83,7 @@ extern int pj_errno; /* global error return code */ #ifdef PROJ_API_H_NOT_INVOKED_AS_PRIMARY_API /* These make the function declarations below conform with classic proj */ typedef PJ *projPJ; /* projPJ is a pointer to PJ */ - typedef projCtx_t *projCtx; /* projCtx is a pointer to projCtx_t */ + typedef struct projCtx_t *projCtx; /* projCtx is a pointer to projCtx_t */ # define projXY XY # define projLP LP # define projXYZ XYZ -- cgit v1.2.3 From d7fda7b431d7af276d40550630390d5e28461140 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Wed, 29 Nov 2017 20:34:08 +0100 Subject: Remove the last remains of the internal selftest --- CMakeLists.txt | 9 --------- appveyor.yml | 4 ++-- nmake.opt | 4 ---- travis/install.sh | 4 ++-- travis/linux_gcc/before_install.sh | 2 +- 5 files changed, 5 insertions(+), 18 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f489f4ed..f6147c3e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -57,15 +57,6 @@ include(Proj4Config) include(Proj4Mac) include(policies) -################################################################################# -# Self-test build config -################################################################################# - -option(SELFTEST "Include self-test in build" OFF) -if(SELFTEST) - add_definitions(-DPJ_SELFTEST) -endif(SELFTEST) - ################################################################################# # threading configuration ################################################################################# diff --git a/appveyor.yml b/appveyor.yml index f8bbd7ef..7de5d23c 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -31,10 +31,10 @@ build_script: - if "%BUILD_TYPE%" == "nmake" if "%platform%" == "x86" call "%VS120COMNTOOLS%\..\..\VC\vcvarsall.bat" - if "%BUILD_TYPE%" == "nmake" if "%platform%" == "x64" "C:\Program Files\Microsoft SDKs\Windows\v7.1\Bin\SetEnv.cmd" /x64 - if "%BUILD_TYPE%" == "nmake" if "%platform%" == "x64" call "%VS120COMNTOOLS%\..\..\VC\vcvarsall.bat" x86_amd64 - - if "%BUILD_TYPE%" == "nmake" nmake /f makefile.vc SELFTEST=1 + - if "%BUILD_TYPE%" == "nmake" nmake /f makefile.vc - if "%BUILD_TYPE%" == "nmake" nmake /f makefile.vc install-all - if "%BUILD_TYPE%" == "nmake" if "%platform%" == "x64" cd src - - if "%BUILD_TYPE%" == "nmake" if "%platform%" == "x64" nmake /f makefile.vc SELFTEST=1 multistresstest.exe + - if "%BUILD_TYPE%" == "nmake" if "%platform%" == "x64" nmake /f makefile.vc multistresstest.exe # Disabled for now as it scales badly # - if "%BUILD_TYPE%" == "nmake" if "%platform%" == "x64" multistresstest.exe - if "%BUILD_TYPE%" == "cmake" if "%platform%" == "x64" SET VS_FULL=%VS_VERSION% Win64 diff --git a/nmake.opt b/nmake.opt index 4e4262fd..1d5b54a7 100644 --- a/nmake.opt +++ b/nmake.opt @@ -28,10 +28,6 @@ OPTFLAGS= /Zi /MDd /Fdproj.pdb !ENDIF !ENDIF -!IFDEF SELFTEST -OPTFLAGS= $(OPTFLAGS) -DPJ_SELFTEST -!ENDIF - # Uncomment the first for linking exes against DLL or second for static EXE_PROJ = proj_i.lib #EXE_PROJ = proj.lib diff --git a/travis/install.sh b/travis/install.sh index 78da24de..8126182b 100755 --- a/travis/install.sh +++ b/travis/install.sh @@ -74,9 +74,9 @@ PROJ_LIB=../nad src/multistresstest cd .. # autoconf build with grids and coverage if [ $TRAVIS_OS_NAME == "osx" ]; then - CFLAGS="-DPJ_SELFTEST --coverage" ./configure; + CFLAGS="--coverage" ./configure; else - CFLAGS="-DPJ_SELFTEST --coverage" LDFLAGS="-lgcov" ./configure; + CFLAGS="--coverage" LDFLAGS="-lgcov" ./configure; fi make -j3 make check diff --git a/travis/linux_gcc/before_install.sh b/travis/linux_gcc/before_install.sh index 0c37643c..87877629 100755 --- a/travis/linux_gcc/before_install.sh +++ b/travis/linux_gcc/before_install.sh @@ -4,7 +4,7 @@ sudo apt-get install -y cppcheck -cppcheck --inline-suppr --template='{file}:{line},{severity},{id},{message}' --enable=all --inconclusive --std=posix -DPJ_SELFTEST=1 src/*.c 2>/tmp/cppcheck.txt +cppcheck --inline-suppr --template='{file}:{line},{severity},{id},{message}' --enable=all --inconclusive --std=posix src/*.c 2>/tmp/cppcheck.txt grep "error," /tmp/cppcheck.txt if [[ $? -eq 0 ]] ; then -- cgit v1.2.3 From a4308d9c529fd913cd8d50813fc913bf80552a38 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Sun, 3 Dec 2017 17:06:02 -0500 Subject: Fixes to geodesic documention: * replace SVG figure by PNG version, * rewrite maths using :math:, * put references in the main reference section, * reformat references for consistency, * put references in alphabetical order, * use USGS URL for Snyder (1987). --- docs/source/geodesic.rst | 187 +++++++++++++++++---------------------------- docs/source/references.rst | 134 +++++++++++++++++++++++++------- 2 files changed, 179 insertions(+), 142 deletions(-) diff --git a/docs/source/geodesic.rst b/docs/source/geodesic.rst index 29327ae8..d54212ca 100644 --- a/docs/source/geodesic.rst +++ b/docs/source/geodesic.rst @@ -10,53 +10,52 @@ Geodesic calculations Introduction ------------ -Consider a ellipsoid of revolution with equatorial radius *a*, polar -semi-axis *b*, and flattening *f* = (*a* − *b*)/*a* . Points on -the surface of the ellipsoid are characterized by their latitude φ -and longitude λ. (Note that latitude here means the +Consider a ellipsoid of revolution with equatorial radius :math:`a`, polar +semi-axis :math:`b`, and flattening :math:`f=(a−b)/a`. Points on +the surface of the ellipsoid are characterized by their latitude :math:`\phi` +and longitude :math:`\lambda`. (Note that latitude here means the *geographical latitude*, the angle between the normal to the ellipsoid and the equatorial plane). The shortest path between two points on the ellipsoid at -(φ\ :sub:`1`, λ\ :sub:`1`) and (φ\ :sub:`2`, λ\ :sub:`2`) +:math:`(\phi_1,\lambda_1)` and :math:`(\phi_2,\lambda_2)` is called the geodesic. Its length is -*s*\ :sub:`12` and the geodesic from point 1 to point 2 has forward -azimuths α\ :sub:`1` and α\ :sub:`2` at the two end -points. In this figure, we have λ\ :sub:`12` = -λ\ :sub:`2` − λ\ :sub:`1`. +:math:`s_{12}` and the geodesic from point 1 to point 2 has forward +azimuths :math:`\alpha_1` and :math:`\alpha_2` at the two end +points. In this figure, we have :math:`\lambda_{12}=\lambda_2-\lambda_1`. .. raw:: html
- Figure from wikipedia + width="320">
A geodesic can be extended indefinitely by requiring that any sufficiently small segment is a shortest path; geodesics are also the straightest curves on the surface. -Solution of geodesic programs +Solution of geodesic problems ----------------------------- Traditionally two geodesic problems are considered: -* the direct problem — given φ\ :sub:`1`, - λ\ :sub:`1`, α\ :sub:`1`, *s*\ :sub:`12`, - determine φ\ :sub:`2`, λ\ :sub:`2`, and - α\ :sub:`2`. +* the direct problem — given :math:`\phi_1`, + :math:`\lambda_1`, :math:`\alpha_1`, :math:`s_{12}`, + determine :math:`\phi_2`, :math:`\lambda_2`, :math:`\alpha_2`. -* the inverse problem — given φ\ :sub:`1`, - λ\ :sub:`1`, φ\ :sub:`2`, λ\ :sub:`2`, - determine *s*\ :sub:`12`, α\ :sub:`1`, and - α\ :sub:`2`. +* the inverse problem — given :math:`\phi_1`, + :math:`\lambda_1`, :math:`\phi_2`, :math:`\lambda_2`, + determine :math:`s_{12}`, :math:`\alpha_1`, + :math:`\alpha_2`. PROJ incorporates `C library for Geodesics `_ from `GeographicLib `_. This library provides routines to solve the direct and inverse geodesic problems. Full double -precision accuracy is maintained provided that −0.02 < *f* < 0.02. Refer +precision accuracy is maintained provided that +:math:`\lvert f\rvert<\frac1{50}`. Refer to the `application programming interface @@ -70,9 +69,9 @@ rest of PROJ: * angles (latitudes, longitudes, and azimuths) are in degrees (instead of in radians); -* the shape of ellipsoid is specified by the flattening *f*; this can - be negative to denote a prolate ellipsoid; setting *f* = 0 corresponds - to a sphere in which case the geodesic becomes a great circle. +* the shape of ellipsoid is specified by the flattening :math:`f`; this can + be negative to denote a prolate ellipsoid; setting :math:`f=0` corresponds + to a sphere, in which case the geodesic becomes a great circle. PROJ also includes a command line tool, :ref:`geod`\ (1), for performing simple geodesic calculations. @@ -82,47 +81,41 @@ Additional properties The routines also calculate several other quantities of interest -* *S*\ :sub:`12` is the area between the geodesic from point 1 to +* :math:`S_{12}` 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 - (φ\ :sub:`1`,λ\ :sub:`1`), (0,λ\ :sub:`1`), - (0,λ\ :sub:`2`), and - (φ\ :sub:`2`,λ\ :sub:`2`). It is given in + :math:`(\phi_1,\lambda_1)`, :math:`(0,\lambda_1)`, + :math:`(0,\lambda_2)`, and + :math:`(\phi_2,\lambda_2)`. It is given in meters\ :sup:`2`. -* *m*\ :sub:`12`, the reduced length of the geodesic is defined such - that if the initial azimuth is perturbed by *d*\ α\ :sub:`1` - (radians) then the second point is displaced by *m*\ :sub:`12` - *d*\ α\ :sub:`1` in the direction perpendicular to the - geodesic. *m*\ :sub:`12` is given in meters. On a curved surface - the reduced length obeys a symmetry relation, *m*\ :sub:`12` + - *m*\ :sub:`21` = 0. On a flat surface, we have *m*\ :sub:`12` = - *s*\ :sub:`12`. -* *M*\ :sub:`12` and *M*\ :sub:`21` are geodesic scales. If two +* :math:`m_{12}`, the reduced length of the geodesic is defined such + that if the initial azimuth is perturbed by :math:`d\alpha_1` + (radians) then the second point is displaced by :math:`m_{12}\,d\alpha_1` + in the direction perpendicular to the + geodesic. :math:`m_{12}` is given in meters. On a curved surface + the reduced length obeys a symmetry relation, :math:`m_{12}+m_{21}=0`. + On a flat surface, we have :math:`m_{12}=s_{12}`. +* :math:`M_{12}` and :math:`M_{21}` 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 *M*\ :sub:`12` *dt* at - point 2. *M*\ :sub:`21` is defined similarly (with the geodesics - being parallel to one another at point 2). *M*\ :sub:`12` and - *M*\ :sub:`21` are dimensionless quantities. On a flat surface, - we have *M*\ :sub:`12` = *M*\ :sub:`21` = 1. -* σ\ :sub:`12` is the arc length on the auxiliary sphere. + :\math`dt`, then they are separated by a distance :math:`M_{12}\,dt` at + point 2. :math:`M_{21}` is defined similarly (with the geodesics + being parallel to one another at point 2). :math:`M_{12}` and + :math:`M_{21}` are dimensionless quantities. On a flat surface, + we have :math:`M_{12}=M_{21}=1`. +* :math:`\sigma_{12}` is the arc length on the auxiliary sphere. This is a construct for converting the problem to one in spherical trigonometry. The spherical arc length from one equator crossing to - the next is always 180°. + the next is always :math:`180^\circ`. If points 1, 2, and 3 lie on a single geodesic, then the following addition rules hold: -* *s*\ :sub:`13` = *s*\ :sub:`12` + *s*\ :sub:`23` -* σ\ :sub:`13` = σ\ :sub:`12` + σ\ :sub:`23` -* *S*\ :sub:`13` = *S*\ :sub:`12` + *S*\ :sub:`23` -* *m*\ :sub:`13` = *m*\ :sub:`12`\ *M*\ :sub:`23` + - *m*\ :sub:`23`\ *M*\ :sub:`21` -* *M*\ :sub:`13` = *M*\ :sub:`12`\ *M*\ :sub:`23` − - (1 − *M*\ :sub:`12`\ *M*\ :sub:`21`) - *m*\ :sub:`23`/*m*\ :sub:`12` -* *M*\ :sub:`31` = *M*\ :sub:`32`\ *M*\ :sub:`21` − - (1 − *M*\ :sub:`23`\ *M*\ :sub:`32`) - *m*\ :sub:`12`/*m*\ :sub:`23` +* :math:`s_{13}=s_{12}+s_{23}`, +* :math:`\sigma_{13}=\sigma_{12}+\sigma_{23}`, +* :math:`S_{13}=S_{12}+S_{23}`, +* :math:`m_{13}=m_{12}M_{23}+m_{23}M_{21}`, +* :math:`M_{13}=M_{12}M_{23}-(1-M_{12}M_{21})m_{23}/m_{12}`, +* :math:`M_{31}=M_{32}M_{21}-(1-M_{23}M_{32})m_{12}/m_{23}`. Multiple shortest geodesics --------------------------- @@ -132,42 +125,41 @@ The shortest distance found by solving the inverse problem is multiple azimuths which yield the same shortest distance. Here is a catalog of those cases: -* φ\ :sub:`1` = −φ\ :sub:`2` (with neither point at - a pole). If α\ :sub:`1` = α\ :sub:`2`, the geodesic +* :math:`\phi_1=-\phi_2` (with neither point at + a pole). If :math:`\alpha_1=\alpha_2`, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting - [α\ :sub:`1`,α\ :sub:`2`] ← [α\ :sub:`2`,α\ :sub:`1`], - [*M*\ :sub:`12`,\ *M*\ :sub:`21`] ← [*M*\ :sub:`21`,\ *M*\ :sub:`12`], - *S*\ :sub:`12` ← −\ *S*\ :sub:`12`. - (This occurs when the longitude difference is near ±180° for oblate - ellipsoids.) -* λ\ :sub:`2` = λ\ :sub:`1` ± 180° (with - neither point at a pole). If α\ :sub:`1` = 0° or - ±180°, the geodesic is unique. Otherwise there are two + :math:`[\alpha_1,\alpha_2]\leftarrow[\alpha_2,\alpha_1]`, + :math:`[M_{12},M_{21}]\leftarrow[M_{21},M_{12}]`, + :math:`S_{12}\leftarrow-S_{12}`. + (This occurs when the longitude difference is near :math:`\pm180^\circ` + for oblate ellipsoids.) +* :math:`\lambda_2=\lambda_1\pm180^\circ` (with + neither point at a pole). If :math:`\alpha_1=0^\circ` or + :math:`\pm180^\circ`, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting - [α\ :sub:`1`,α\ :sub:`2`] ← [−α\ :sub:`1`,−α\ :sub:`2`], - *S*\ :sub:`12` ← −\ *S*\ :sub:`12`. (This occurs when - φ\ :sub:`2` is near −φ\ :sub:`1` for prolate + :math:`[\alpha_1,\alpha_2]\leftarrow[-\alpha_1,-\alpha_2]`, + :math:`S_{12}\leftarrow-S_{12}`. (This occurs when + :math:`\phi_2` is near :math:`-\phi_1` for prolate ellipsoids.) * Points 1 and 2 at opposite poles. There are infinitely many geodesics which can be generated by setting - [α\ :sub:`1`,α\ :sub:`2`] ← - [α\ :sub:`1`,α\ :sub:`2`] + [δ,−δ], for arbitrary δ. + :math:`[\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,-\delta]`, + for arbitrary :math:`\delta`. (For spheres, this prescription applies when points 1 and 2 are antipodal.) -* *s*\ :sub:`12` = 0 (coincident points). There are infinitely many +* :math:`s_{12}=0` (coincident points). There are infinitely many geodesics which can be generated by setting - [α\ :sub:`1`,α\ :sub:`2`]_← - [α\ :sub:`1`,α\ :sub:`2`]_+ [δ,δ], for - arbitrary δ. + :math:`[\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,\delta]`, + for arbitrary :math:`\delta`. Background ---------- -The algorithms implemented by this package are given in Karney (2013) -and are based on Bessel (1825) and Helmert (1880); the algorithm for -areas is based on Danielsen (1989). These improve on the work of -Vincenty (1975) in the following respects: +The algorithms implemented by this package are given in [Karney2013]_ +and are based on [Bessel1825]_ and [Helmert1880]_; the algorithm for +areas is based on [Danielsen1989]_. These improve on the work of +[Vincenty1975]_ in the following respects: * The results are accurate to round-off for terrestrial ellipsoids (the error in the distance is less then 15 nanometers, compared to 0.1 mm @@ -178,40 +170,5 @@ Vincenty (1975) in the following respects: geodesic. This allows, for example, the area of a geodesic polygon to be computed. -References ----------- - -* F. W. Bessel, - `The calculation of longitude and latitude from geodesic measurements (1825) - `_, - Astron. Nachr. **331**\ (8), 852–861 (2010), - translated by C. F. F. Karney and R. E. Deakin. -* F. R. Helmert, - `Mathematical and Physical Theories of Higher Geodesy, Vol 1 - `_, - (Teubner, Leipzig, 1880), Chaps. 5–7. -* T. Vincenty, - `Direct and inverse solutions of geodesics on the ellipsoid with - application of nested equations - `_, - Survey Review **23**\ (176), 88–93 (1975). -* J. Danielsen, - `The area under the geodesic - `_, - Survey Review **30**\ (232), 61–66 (1989). -* C. F. F. Karney, - `Algorithms for geodesics - `_, - J. Geodesy **87**\ (1) 43–55 (2013); - `addenda `_. -* C. F. F. Karney, - `Geodesics on an ellipsoid of revolution - `_, - Feb. 2011; - `errata - `_. -* `A geodesic bibliography - `_. -* The wikipedia page, - `Geodesics on an ellipsoid - `_. +Additional background material is provided in [GeodesicBib]_, +[GeodesicWiki]_, and [Karney2011]_. diff --git a/docs/source/references.rst b/docs/source/references.rst index 31100ded..e80e8caf 100644 --- a/docs/source/references.rst +++ b/docs/source/references.rst @@ -5,30 +5,110 @@ References ================================================================================ -.. [Evenden1995] Evenden, G. I., 1995, `Cartograpic Projection Procedures for the UNIX Environment - A User's Manual `_ - -.. [Evenden2005] Evenden, G. I., 2005, `libproj4: A Comprehensive Library of Cartographic Projection Functions (Preliminary Draft) `_ - -.. [Steers1970] Steers, J.A., 1970, An introduction to the study of map projections (15th ed.): London, Univ. London Press, p. 229 - -.. [Snyder1987] Snyder. John P. 1987. `Map Projections - A Working Manual `_. US. Geological Survey professional paper; 1395. - -.. [Snyder1993] Snyder, 1993, Flattening the Earth, Chicago and London, The university of Chicago press - -.. [EberHewitt1979] Eber, L.E., and R.P. Hewitt. 1979. `Conversion algorithms for the CALCOFI station grid `__. California Cooperative Oceanic Fisheries Investigations Reports 20:135-137. - -.. [WeberMoore2013] Weber, E.D., and T.J. Moore. 2013. `Corrected Conversion Algorithms For The Calcofi Station Grid And Their Implementation In Several Computer Languages `__. California Cooperative Oceanic Fisheries Investigations Reports 54. - -.. [CalabrettaGreisen2002] M. Calabretta and E. Greisen, 2002, "Representations of celestial coordinates in FITS". Astronomy & Astrophysics 395, 3, 1077–1122. - -.. [ChanONeil1975] F. Chan and E.M.O'Neill, 1975, "Feasibility Study of a Quadrilateralized Spherical Cube Earth Data Base". Tech. Rep. EPRF 2-75 (CSC), Environmental Prediction Research Facility. - -.. [ONeilLaubscher1976] E.M. O'Neill and R.E. Laubscher, 1976, "Extended Studies of a Quadrilateralized Spherical Cube Earth Data Base". Tech. Rep. NEPRF 3-76 (CSC), Naval Environmental Prediction Research Facility. - -.. [LambersKolb2012] M. Lambers and A. Kolb, 2012, "Ellipsoidal Cube Maps for Accurate Rendering of Planetary-Scale Terrain Data", Proc. Pacfic Graphics (Short Papers). - -.. [Zajac1978] A. Zajac, 1978, "Atlas of distribution of vascular plants in Poland (ATPOL)". Taxon 27(5/6), 481–484. - -.. [Komsta2016] L. Komsta, 2016, `ATPOL geobotanical grid revisited – a proposal of coordinate conversion algorithms `__. Annales UMCS Sectio E Agricultura 71(1), 31-37. - -.. [Verey2017] M. Verey, 2017, `Theoretical analysis and practical consequences of adopting an ATPOL grid model as a conical projection, defining the conversion of plane coordinates to the WGS - 84 ellipsoid `__. Fragmenta Floristica et Geobotanica Polonica (preprint submitted). +.. [Bessel1825] F. W. Bessel, 1825, + `The calculation of longitude and latitude from geodesic measurements + `_, + Astron. Nachr. **331**\ (8), 852–861 (2010), + translated by C. F. F. Karney and R. E. Deakin. + +.. [CalabrettaGreisen2002] M. Calabretta and E. Greisen, 2002, + "Representations of celestial coordinates in FITS". + Astronomy & Astrophysics 395, 3, 1077–1122. + +.. [ChanONeil1975] F. Chan and E. M. O'Neill, 1975, + "Feasibility Study of a Quadrilateralized Spherical Cube Earth Data Base", + Tech. Rep. EPRF 2-75 (CSC), Environmental Prediction Research Facility. + +.. [Danielsen1989] J. Danielsen, 1989, + `The area under the geodesic + `_, + Survey Review **30**\ (232), 61–66. + +.. [EberHewitt1979] L. E. Eber and R.P. Hewitt, 1979, + `Conversion algorithms for the CALCOFI station grid + `_, + California Cooperative Oceanic Fisheries Investigations Reports 20:135-137. + +.. [Evenden1995] G. I. Evenden, 1995, + `Cartograpic Projection Procedures for the UNIX Environment - + A User's Manual + `_. + +.. [Evenden2005] G. I. Evenden, 2005, + `libproj4: A Comprehensive Library of Cartographic Projection Functions + (Preliminary Draft) + `_. + +.. [GeodesicBib] `A geodesic bibliography + `_. + +.. [GeodesicWiki] The wikipedia page, + `Geodesics on an ellipsoid + `_. + +.. [Helmert1880] F. R. Helmert, 1880, + `Mathematical and Physical Theories of Higher Geodesy, Vol 1 + `_, + (Teubner, Leipzig), Chaps. 5–7. + +.. [Karney2011] C. F. F. Karney, 2011, + `Geodesics on an ellipsoid of revolution + `_; + `errata + `_. + +.. [Karney2013] C. F. F. Karney, 2013, + `Algorithms for geodesics + `_, + J. Geodesy **87**\ (1) 43–55; + `addenda `_. + +.. [Komsta2016] L. Komsta, 2016, + `ATPOL geobotanical grid revisited - a proposal of coordinate conversion + algorithms + `_, + Annales UMCS Sectio E Agricultura 71(1), 31–37. + +.. [LambersKolb2012] M. Lambers and A. Kolb, 2012, + "Ellipsoidal Cube Maps for Accurate Rendering of Planetary-Scale + Terrain Data", Proc. Pacfic Graphics (Short Papers). + +.. [ONeilLaubscher1976] E. M. O'Neill and R. E. Laubscher, 1976, + "Extended Studies of a Quadrilateralized Spherical Cube Earth Data Base", + Tech. Rep. NEPRF 3-76 (CSC), + Naval Environmental Prediction Research Facility. + +.. [Snyder1987] J. P. Snyder, 1987, + `Map Projections - A Working Manual + `_. + U.S. Geological Survey professional paper 1395. + +.. [Snyder1993] J. P. Snyder, 1993, + Flattening the Earth, Chicago and London, The University of Chicago press. + +.. [Steers1970] J. A. Steers, 1970, + An introduction to the study of map projections (15th ed.), + London Univ. Press, p. 229. + +.. [Verey2017] M. Verey, 2017, + `Theoretical analysis and practical consequences of adopting an ATPOL + grid model as a conical projection, defining the conversion of plane + coordinates to the WGS-84 ellipsoid + `_, + Fragmenta Floristica et Geobotanica Polonica (preprint submitted). + +.. [Vincenty1975] T. Vincenty, 1975, + `Direct and inverse solutions of geodesics on the ellipsoid with + application of nested equations + `_, + Survey Review **23**\ (176), 88–93. + +.. [WeberMoore2013] E. D. Weber and T.J. Moore, 2013, + `Corrected Conversion Algorithms For The Calcofi Station Grid And Their + Implementation In Several Computer Languages + `_, + California Cooperative Oceanic Fisheries Investigations Reports 54. + +.. [Zajac1978] A. Zajac, 1978, + "Atlas of distribution of vascular plants in Poland (ATPOL)", + Taxon 27(5/6), 481–484. -- cgit v1.2.3