Merge remote-tracking branch 'osgeo/master' into docs-release-4.10.0

27 files changed, 459 insertions, 397 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index f489f4ed..f6147c3e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -58,15 +58,6 @@ 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
 #################################################################################
 set(CMAKE_THREAD_PREFER_PTHREAD TRUE)
diff --git a/appveyor.yml b/appveyor.yml
index 43669a40..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
@@ -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/docs/source/faq.rst b/docs/source/faq.rst
index 48bcfafd..8afdb35d 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..d54212ca 100644
--- a/docs/source/geodesic.rst
+++ b/docs/source/geodesic.rst
@@ -1,219 +1,174 @@
 .. _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 <https://en.wikipedia.org/wiki/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) <http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf>`__
-
-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 <http://geographiclib.sourceforge.net/scripts/geod-calc.html>`__,
-   a text interface to geodesic calculations;
-* `geod-google <http://geographiclib.sourceforge.net/scripts/geod-google.html>`__,
-   a tool for drawing geodesics on Google Maps.
-
-Implementations in `Python <http://pypi.python.org/pypi/geographiclib>`__,
-`Matlab <http://www.mathworks.com/matlabcentral/fileexchange/39108>`__,
-`C <http://geographiclib.sourceforge.net/html/C/>`__,
-`Fortran <http://geographiclib.sourceforge.net/html/Fortran/>`__ , and
-`Java <http://geographiclib.sourceforge.net/html/java/>`__ are also available.
-
-The algorithms are described in
- * C. F. F. Karney, `Algorithms for gedesics <http://dx.doi.org/10.1007/s00190-012-0578-z>`__,
-   J. Geodesy '''87'''(1), 43-55 (2013),
-   DOI: `10.1007/s00190-012-0578-z <http://dx.doi.org/10.1007/s00190-012-0578-z>`__; `geo-addenda.html <http://geographiclib.sf.net/geod-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 <https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid#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 :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
+:math:`(\phi_1,\lambda_1)` and :math:`(\phi_2,\lambda_2)`
+is called the geodesic.  Its length is
+: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
+
+       <center>
+         <img src="https://upload.wikimedia.org/wikipedia/commons/thumb/c/cb/Geodesic_problem_on_an_ellipsoid.svg/320px-Geodesic_problem_on_an_ellipsoid.svg.png"
+              alt="Figure from wikipedia"
+              width="320">
+       </center>
+
+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 problems
+-----------------------------
+
+Traditionally two geodesic problems are considered:
+
+* 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  :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
+<https://geographiclib.sourceforge.io/1.49/C/>`_ from `GeographicLib
+<https://geographiclib.sourceforge.io>`_.  This library provides
+routines to solve the direct and inverse geodesic problems.  Full double
+precision accuracy is maintained provided that
+:math:`\lvert f\rvert<\frac1{50}`.  Refer
+to the
+
+    `application programming interface
+    <https://geographiclib.sourceforge.io/1.49/C/geodesic_8h.html>`_
+
+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 :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.
+
+Additional properties
+---------------------
+
+The routines also calculate several other quantities of interest
+
+* :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
+  :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`.
+* :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
+  :\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 :math:`180^\circ`.
+
+If points 1, 2, and 3 lie on a single geodesic, then the following
+addition rules hold:
+
+* :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
+---------------------------
+
+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:
+
+* :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
+  :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
+  :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
+  :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.)
+* :math:`s_{12}=0` (coincident points).  There are infinitely many
+  geodesics which can be generated by setting
+  :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 [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
+  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.
+
+Additional background material is provided in [GeodesicBib]_,
+[GeodesicWiki]_, and [Karney2011]_.
diff --git a/docs/source/references.rst b/docs/source/references.rst
index e2dece5c..ad5f16bf 100644
--- a/docs/source/references.rst
+++ b/docs/source/references.rst
@@ -4,37 +4,125 @@
 References
 ================================================================================
 
-.. [AltamimiEtAl2002] Altamimi, Z., P. Sillard, and C. Boucher (2002), ITRF2000: A new release of the International Terrestrial Reference Frame for earth science applications, J. Geophys. Res., 107(B10), 2214, `doi:10.1029/2001JB000561 <http://www.agu.org/pubs/crossref/2002/2001JB000561.shtml>`__.
-
-.. [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.
-
-.. [Deakin2004] R.E. Deakin, 2004, `The Standard and Abridged Molodensky Coordinate Transformation Formulae <http://www.mygeodesy.id.au/documents/Molodensky%20V2.pdf>`__.
-
-.. [EberHewitt1979] Eber, L.E., and R.P. Hewitt. 1979. `Conversion algorithms for the CALCOFI station grid <http://www.calcofi.org/publications/calcofireports/v20/Vol_20_Eber___Hewitt.pdf>`__. California Cooperative Oceanic Fisheries Investigations Reports 20:135-137.
-
-.. [Evenden1995] Evenden, G. I., 1995, `Cartograpic Projection Procedures for the UNIX Environment - A User's Manual <https://github.com/OSGeo/proj.4/blob/master/docs/old/proj_4_3_12.pdf>`_
-
-.. [Evenden2005] Evenden, G. I., 2005, `libproj4: A Comprehensive Library of Cartographic Projection Functions (Preliminary Draft) <https://github.com/OSGeo/proj.4/blob/master/docs/old/libproj.pdf>`_
-
-.. [EversKnudsen2017] K. Evers and T. Knudsen, 2017, `Transformation pipelines for PROJ.4 <http://www.fig.net/resources/proceedings/fig_proceedings/fig2017/papers/iss6b/ISS6B_evers_knudsen_9156.pdf>`__, FIG Working Week 2017 Proceedings.
-
-.. [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.
-
-.. [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 <https://github.com/OSGeo/proj.4/blob/master/docs/old/USGS-Snyder-Map-Projections-A-Working-Manual-1987.pdf>`_. US. Geological Survey professional paper; 1395.
-
-.. [Snyder1993] Snyder, 1993, Flattening the Earth, Chicago and London, The university of Chicago press
-
-.. [WeberMoore2013] Weber, E.D., and T.J. Moore. 2013. `Corrected Conversion Algorithms For The Calcofi Station Grid And Their Implementation In Several Computer Languages <http://calcofi.org/publications/calcofireports/v54/Vol_54_Weber.pdf>`__. 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.
-
-.. [Komsta2016] L. Komsta, 2016, `ATPOL geobotanical grid revisited – a proposal of coordinate conversion algorithms <http://wydawnictwo.up.lublin.pl/annales/Agricultura/2016/1/03.pdf>`__. 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 <http://www.botany.pl/atpol/Siatka%20ATPOL%20w%20analitycznym%20ujeciu.pdf>`__. Fragmenta Floristica et Geobotanica Polonica (preprint submitted).
+.. [AltamimiEtAl2002] Altamimi, Z., P. Sillard, and C. Boucher (2002),
+    ITRF2000: A new release of the International Terrestrial Reference Frame for earth science applications,
+    J. Geophys. Res., 107(B10), 2214,
+    `doi:10.1029/2001JB000561 <http://www.agu.org/pubs/crossref/2002/2001JB000561.shtml>`__.
+
+
+.. [Bessel1825] F. W. Bessel, 1825,
+   `The calculation of longitude and latitude from geodesic measurements
+   <https://arxiv.org/abs/0908.1824>`_,
+   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
+   <https://doi.org/10.1179/003962689791474267>`_,
+   Survey Review **30**\ (232), 61–66.
+
+.. [Deakin2004] R.E. Deakin, 2004,
+    `The Standard and Abridged Molodensky Coordinate Transformation Formulae
+    <http://www.mygeodesy.id.au/documents/Molodensky%20V2.pdf>`__.
+
+.. [EberHewitt1979] L. E. Eber and R.P. Hewitt, 1979,
+   `Conversion algorithms for the CALCOFI station grid
+   <http://www.calcofi.org/publications/calcofireports/v20/Vol_20_Eber___Hewitt.pdf>`_,
+   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
+   <https://github.com/OSGeo/proj.4/blob/master/docs/old/proj_4_3_12.pdf>`_.
+
+.. [Evenden2005] G. I. Evenden, 2005,
+   `libproj4: A Comprehensive Library of Cartographic Projection Functions
+   (Preliminary Draft)
+   <https://github.com/OSGeo/proj.4/blob/master/docs/old/libproj.pdf>`_.
+
+.. [EversKnudsen2017] K. Evers and T. Knudsen, 2017,
+    `Transformation pipelines for PROJ.4
+    <http://www.fig.net/resources/proceedings/fig_proceedings/fig2017/papers/iss6b/ISS6B_evers_knudsen_9156.pdf>`__,
+    FIG Working Week 2017 Proceedings.
+
+.. [GeodesicBib] `A geodesic bibliography
+   <https://geographiclib.sourceforge.io/geodesic-papers/biblio.html>`_.
+
+.. [GeodesicWiki] The wikipedia page,
+   `Geodesics on an ellipsoid
+   <https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid>`_.
+
+.. [Helmert1880] F. R. Helmert, 1880,
+   `Mathematical and Physical Theories of Higher Geodesy, Vol 1
+   <https://doi.org/10.5281/zenodo.32050>`_,
+   (Teubner, Leipzig), Chaps. 5–7.
+
+.. [Karney2011] C. F. F. Karney, 2011,
+   `Geodesics on an ellipsoid of revolution
+   <https://arxiv.org/abs/1102.1215v1>`_;
+   `errata
+   <https://geographiclib.sourceforge.io/geod-addenda.html#geod-errata>`_.
+
+.. [Karney2013] C. F. F. Karney, 2013,
+   `Algorithms for geodesics
+   <https://doi.org/10.1007/s00190-012-0578-z>`_,
+   J. Geodesy **87**\ (1) 43–55;
+   `addenda <https://geographiclib.sourceforge.io/geod-addenda.html>`_.
+
+.. [Komsta2016] L. Komsta, 2016,
+   `ATPOL geobotanical grid revisited - a proposal of coordinate conversion
+   algorithms
+   <http://wydawnictwo.up.lublin.pl/annales/Agricultura/2016/1/03.pdf>`_,
+   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
+   <https://pubs.er.usgs.gov/publication/pp1395>`_.
+   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
+   <http://www.botany.pl/atpol/Siatka%20ATPOL%20w%20analitycznym%20ujeciu.pdf>`_,
+   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
+   <http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf>`_,
+   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
+   <http://calcofi.org/publications/calcofireports/v54/Vol_54_Weber.pdf>`_,
+   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.
diff --git a/docs/source/usage/apps/geod.rst b/docs/source/usage/apps/geod.rst
index 7765f36a..1e009283 100644
--- a/docs/source/usage/apps/geod.rst
+++ b/docs/source/usage/apps/geod.rst
@@ -160,11 +160,9 @@ which gives:
 Further reading
 ***************
 
-#. `GeographicLib <http://geographiclib.sf.net>`_
-
-#. `C. F. F. Karney, Algorithms for Geodesics, J. Geodesy 87, 43-55 (2013) <http://dx.doi.org/10.1007/s00190-012-0578-z>`_.
-   `Addendum <http://geographiclib.sf.net/geod-addenda.html>`_
-
-#. `The online geodesic bibliography <http://geographiclib.sf.net/geodesic-papers/biblio.html>`_
+#. `GeographicLib <https://geographiclib.sourceforge.io>`_.
 
+#. C. F. F. Karney, `Algorithms for Geodesics <https://doi.org/10.1007/s00190-012-0578-z>`_, J. Geodesy **87**\ (1), 43–55 (2013);
+   `addenda <https://geographiclib.sourceforge.io/geod-addenda.html>`_.
 
+#. `A geodesic bibliography <https://geographiclib.sourceforge.io/geodesic-papers/biblio.html>`_.
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
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/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 <assert.h>
 #include <stddef.h>
 #include <math.h>
 #include <errno.h>
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index d0fd8222..34bb7a68 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -43,7 +43,6 @@ Last update: 2017-05-15
 ***********************************************************************/
 
 #define PJ_LIB__
-#include <assert.h>
 #include <stddef.h>
 #include <errno.h>
 #include "proj_internal.h"
@@ -572,7 +571,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/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 <assert.h>
 #include <stddef.h>
 #include <math.h>
 #include <errno.h>
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 <assert.h>
 #include <stddef.h>
 #include <errno.h>
 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 <stdio.h>
 #include <stdlib.h>
-#include <assert.h>
 
 #define PJ_LIB__
-#include <projects.h>
+#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 <errno.h>
 #include <string.h>
 
-#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 <assert.h>
-#endif /* _WIN32_WCE */
-
 /************************************************************************/
 /*                             swap_words()                             */
 /*                                                                      */
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_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 <projects.h>
+#include "proj_internal.h"
+#include "projects.h"
 #include <string.h>
 #include <math.h>
 #include <errno.h>
 
-#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 <assert.h>
-#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 <string.h>
 #include <math.h>
 
-#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 <assert.h>
-#endif /* _WIN32_WCE */
-
 static PJ_GRIDINFO *grid_list = NULL;
 #define PJ_MAX_PATH_LENGTH 1024
 
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_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/src/proj_api.h b/src/proj_api.h
index b57d9f38..597a2589 100644
--- a/src/proj_api.h
+++ b/src/proj_api.h
@@ -69,22 +69,21 @@ extern "C" {
 
 
 
-#ifndef PROJ_H
 extern char const pj_release[]; /* global release id string */
 extern int pj_errno;    /* global error return code */
 
-/* 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 */
-    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
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 <limits.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -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 */
diff --git a/test/gie/GDA.gie b/test/gie/GDA.gie
new file mode 100644
index 00000000..3fe0f3e3
--- /dev/null
+++ b/test/gie/GDA.gie
@@ -0,0 +1,70 @@
+-----------------------------------------------------------------------------------
+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
+-----------------------------------------------------------------------------------
+
+<gie>
+-----------------------------------------------------------------------------------
+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;   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
+
+ # Alice Springs ITRF2014@2018.0
+accept     -4052052.6588 4212835.9938 -2545104.6946 2018.0
+
+ # Alice Springs GDA2020
+expect     -4052052.7373 4212835.9835 -2545104.5867
+
+</gie>
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
diff --git a/travis/install.sh b/travis/install.sh
index f71108c6..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
@@ -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
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