22 files changed, 864 insertions, 29 deletions
diff --git a/docs/plot/plot.py b/docs/plot/plot.py
index 139281bf..a4f330f6 100644
--- a/docs/plot/plot.py
+++ b/docs/plot/plot.py
@@ -157,8 +157,14 @@ def project(coordinates, proj_string, in_radians=False):
     args = [PROJ, '-b']
     args.extend(proj_string.split(' '))
 
-    proc = subprocess.Popen(args, stdin=subprocess.PIPE, stdout=subprocess.PIPE,
-                            env={'PROJ_LIB': os.path.abspath(PROJ_LIB)})
+    try:
+        proc = subprocess.Popen(args, stdin=subprocess.PIPE, stdout=subprocess.PIPE,
+                                env={'PROJ_LIB': os.path.abspath(PROJ_LIB)})
+    except FileNotFoundError:
+        print("'proj' binary not found, please update PROJ constant in plot.py "
+              "to point to your local 'proj' binary")
+        exit(1)
+
     stdout, _ = proc.communicate(coordinates.tobytes(order='C'))
 
     out = np.frombuffer(stdout, dtype=np.double)
@@ -304,7 +310,7 @@ def plotproj(plotdef, data, outdir):
 
     # Make sure the plot is not stretched
     axes.set_aspect('equal')
-    
+
     if not os.path.exists(outdir):
         os.makedirs(outdir)
     plt.savefig(outdir + '/' + plotdef['filename'],
diff --git a/docs/plot/plotdefs.json b/docs/plot/plotdefs.json
index c494d6a7..daaf2ef8 100644
--- a/docs/plot/plotdefs.json
+++ b/docs/plot/plotdefs.json
@@ -88,6 +88,18 @@
         "type": "poly"
     },
     {
+        "filename": "bertin1953.png",
+        "latmax": 90,
+        "latmin": -90,
+        "latmin": -60,
+        "lonmax": 180,
+        "lonmin": -180,
+        "name": "bertin1953",
+        "projstring": "+proj=bertin1953",
+        "res": "low",
+        "type": "poly"
+    },
+    {
         "filename": "bipc.png",
         "latmax": 90,
         "latmin": -90,
@@ -1243,6 +1255,17 @@
         "type": "poly"
     },
     {
+        "filename": "tobmerc.png",
+        "latmax": 80,
+        "latmin": -80,
+        "lonmax": 180,
+        "lonmin": -180,
+        "name": "tobmerc",
+        "projstring": "+proj=tobmerc",
+        "res": "low",
+        "type": "poly"
+    },
+    {
         "filename": "tpeqd.png",
         "latmax": 90,
         "latmin": -90,
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 8892dffc..434e4d73 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -85,7 +85,7 @@ language = None
 
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
-exclude_patterns = []
+exclude_patterns = ['operations/options/*.rst']
 
 # The reST default role (used for this markup: `text`) to use for all
 # documents.
diff --git a/docs/source/operations/projections/bertin1953.rst b/docs/source/operations/projections/bertin1953.rst
new file mode 100644
index 00000000..8be521f3
--- /dev/null
+++ b/docs/source/operations/projections/bertin1953.rst
@@ -0,0 +1,64 @@
+.. _bertin1953:
+
+********************************************************************************
+Bertin 1953
+********************************************************************************
+
+.. versionadded:: 6.0.0
+
++---------------------+-------------------------------+
+| **Classification**  | Miscellaneous                 |
++---------------------+-------------------------------+
+| **Available forms** | Forward, spherical projection |
++---------------------+-------------------------------+
+| **Defined area**    | Global                        |
++---------------------+-------------------------------+
+| **Alias**           | bertin1953                    |
++---------------------+-------------------------------+
+| **Domain**          | 2D                            |
++---------------------+-------------------------------+
+| **Input type**      | Geodetic coordinates          |
++---------------------+-------------------------------+
+| **Output type**     | Projected coordinates         |
++---------------------+-------------------------------+
+
+
+
+.. figure:: ./images/bertin1953.png
+   :width: 500 px
+   :align: center
+   :alt:   Bertin 1953
+
+   proj-string: ``+proj=bertin1953``
+
+
+The Bertin 1953 projection is intended for making world maps. Created by Jacques Bertin
+in 1953, this projection was the go-to choice of the French cartographic school when they
+wished to represent phenomena on a global scale. The projection was formulated in 2017
+by Philippe Rivière for visionscarto.net.
+
+
+Usage
+###############################################################################
+
+The Bertin 1953 projection has no special options. Its rotation parameters
+are fixed. Here is an example of a forward projection with scale 1:
+
+    $ echo 122 47 | src/proj +proj=bertin1953 +R=1
+    0.72	0.73
+
+Parameters
+################################################################################
+
+.. note:: All parameters for the projection are optional.
+
+.. include:: ../options/R.rst
+
+.. include:: ../options/x_0.rst
+
+.. include:: ../options/y_0.rst
+
+Further reading
+################################################################################
+
+#. Philippe Rivière (2017). `Bertin Projection (1953) <https://visionscarto.net/bertin-projection-1953>`, Visionscarto.net.
diff --git a/docs/source/operations/projections/images/bertin1953.png b/docs/source/operations/projections/images/bertin1953.png
new file mode 100644
index 00000000..9b272033
--- /dev/null
+++ b/docs/source/operations/projections/images/bertin1953.png
diff --git a/docs/source/operations/projections/images/tobmerc.png b/docs/source/operations/projections/images/tobmerc.png
new file mode 100644
index 00000000..291bdac9
--- /dev/null
+++ b/docs/source/operations/projections/images/tobmerc.png
diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst
index a48eabe6..572deb1a 100644
--- a/docs/source/operations/projections/index.rst
+++ b/docs/source/operations/projections/index.rst
@@ -20,6 +20,7 @@ Projections map the spherical 3D space to a flat 2D space.
    apian
    august
    bacon
+   bertin1953
    bipc
    boggs
    bonne
@@ -124,6 +125,7 @@ Projections map the spherical 3D space to a flat 2D space.
    tcea
    tissot
    tmerc
+   tobmerc
    tpeqd
    tpers
    ups
diff --git a/docs/source/operations/projections/lcc.rst b/docs/source/operations/projections/lcc.rst
index 6db2a73c..3b2032d6 100644
--- a/docs/source/operations/projections/lcc.rst
+++ b/docs/source/operations/projections/lcc.rst
@@ -4,6 +4,39 @@
 Lambert Conformal Conic
 ********************************************************************************
 
+A Lambert Conformal Conic projection (LCC) is a conic map projection
+used for aeronautical charts, portions of the State Plane Coordinate
+System, and many national and regional mapping systems. It is one of
+seven projections introduced by Johann Heinrich Lambert in 1772.
+
+It has several different forms: with one and two standard parallels
+(referred to as 1SP and 2SP in EPSG guidance notes). Additionally we
+provide "2SP Michigan" form which is very similar to normal 2SP, but
+with a scaling factor on the ellipsoid (given as `k_0` parameter).
+It is implemented as per EPSG Guidance Note 7-2 (version 54, August
+2018, page 25). It is used in a few systems in the EPSG database which
+justifies adding this otherwise non-standard projection.
+
++---------------------+----------------------------------------------------------+
+| **Classification**  | Conformal conic                                          |
++---------------------+----------------------------------------------------------+
+| **Available forms** | Forward and inverse, spherical and elliptical projection.|
+|                     | One or two standard parallels (1SP and 2SP).             |
+|                     | "LCC 2SP Michigan" form can be used by setting           |
+|                     | `+k_0` parameter to specify elliposid scale.             |
++---------------------+----------------------------------------------------------+
+| **Defined area**    | Best for regions predominantly east–west in extent and   |
+|                     | located in the middle north or south latitudes.          |
++---------------------+----------------------------------------------------------+
+| **Alias**           | lcc                                                      |
++---------------------+----------------------------------------------------------+
+| **Domain**          | 2D                                                       |
++---------------------+----------------------------------------------------------+
+| **Input type**      | Geodetic coordinates                                     |
++---------------------+----------------------------------------------------------+
+| **Output type**     | Projected coordinates                                    |
++---------------------+----------------------------------------------------------+
+
 .. figure:: ./images/lcc.png
    :width: 500 px
    :align: center
@@ -31,3 +64,21 @@ Parameters
 .. include:: ../options/x_0.rst
 
 .. include:: ../options/y_0.rst
+
+.. option:: +k_0=<value>
+
+    This parameter can represent two different values depending on the
+    form of the projection. In LCC 1SP it determines the scale factor
+    at natural origin. In LCC 2SP Michigan it determines the ellipsoid
+    scale factor.
+
+    *Defaults to 1.0.*
+
+Further reading
+###############
+
+#. `Wikipedia <https://en.wikipedia.org/wiki/Lambert_conformal_conic_projection>`_
+#. `Wolfram Mathworld <http://mathworld.wolfram.com/LambertConformalConicProjection.html>`_
+#. `John P. Snyder "Map projections: A working manual" (pp. 104-110) <https://pubs.er.usgs.gov/publication/pp1395>`_
+#. `ArcGIS documentation on "Lambert Conformal Conic" <http://desktop.arcgis.com/en/arcmap/10.3/guide-books/map-projections/lambert-conformal-conic.htm>`_
+#. `EPSG Guidance Note 7-2 (version 54, August 2018, page 25) <http://www.epsg.org/Guidancenotes.aspx>`_
diff --git a/docs/source/operations/projections/tobmerc.rst b/docs/source/operations/projections/tobmerc.rst
new file mode 100644
index 00000000..d9a3eb05
--- /dev/null
+++ b/docs/source/operations/projections/tobmerc.rst
@@ -0,0 +1,100 @@
+.. _tobmerc:
+
+********************************************************************************
+Tobler-Mercator
+********************************************************************************
+
+.. versionadded:: 6.0.0
+
+Equal area cylindrical projection with the same latitudinal spacing as
+Mercator projection.
+
++---------------------+----------------------------------------------------------+
+| **Classification**  | Cylindrical equal area                                   |
++---------------------+----------------------------------------------------------+
+| **Available forms** | Forward and inverse, spherical only                      |
++---------------------+----------------------------------------------------------+
+| **Defined area**    | Global, conventionally truncated at about 80 degrees     |
+|                     | north and south                                          |
++---------------------+----------------------------------------------------------+
+| **Alias**           | tobmerc                                                  |
++---------------------+----------------------------------------------------------+
+| **Domain**          | 2D                                                       |
++---------------------+----------------------------------------------------------+
+| **Input type**      | Geodetic coordinates                                     |
++---------------------+----------------------------------------------------------+
+| **Output type**     | Projected coordinates                                    |
++---------------------+----------------------------------------------------------+
+
+.. figure:: ./images/tobmerc.png
+   :width: 500 px
+   :align: center
+   :alt:   Tobler-Mercator
+
+   proj-string: ``+proj=tobmerc``
+
+Usage
+#####
+
+The inappropriate use of the Mercator projection has declined but still
+occasionally occurs. One method of contrasting the Mercator projection is to
+present an alternative in the form of an equal area projection. The map
+projection derived here is thus not simply a pretty Christmas tree ornament:
+it is instead a complement to Mercator's conformal navigation anamorphose
+and can be displayed as an alternative. The equations for the new map
+projection preserve the latitudinal stretching of the Mercator while
+adjusting the longitudinal spacing. This allows placement of the new map
+adjacent to that of Mercator. The surface area, while drastically warped,
+maintains the correct magnitude.
+
+Parameters
+################################################################################
+
+.. note:: All parameters for the projection are optional.
+
+.. include:: ../options/k_0.rst
+
+.. include:: ../options/lon_0.rst
+
+.. include:: ../options/x_0.rst
+
+.. include:: ../options/y_0.rst
+
+.. include:: ../options/R.rst
+
+Mathematical definition
+#######################
+
+The formulas describing the Tobler-Mercator are taken from Waldo Tobler's
+article :cite:`Tobler2017`
+
+Spherical form
+**************
+For the spherical form of the projection we introduce the scaling factor:
+
+.. math::
+
+    k_0 = \cos^2 \phi_{ts}
+
+Forward projection
+==================
+
+.. math::
+
+    x = k_0 \lambda
+
+.. math::
+
+    y = k_0 \ln \left[ \tan \left(\frac{\pi}{4} + \frac{\phi}{2} \right) \right]
+
+
+Inverse projection
+==================
+
+.. math::
+
+    \lambda = \frac{x}{k_0}
+
+.. math::
+
+    \phi = \frac{\pi}{2} - 2 \arctan \left[ e^{-y/k_0} \right]
diff --git a/docs/source/operations/transformations/horner.rst b/docs/source/operations/transformations/horner.rst
new file mode 100644
index 00000000..dd9b9cc9
--- /dev/null
+++ b/docs/source/operations/transformations/horner.rst
@@ -0,0 +1,213 @@
+.. _horner:
+
+================================================================================
+Horner polynomial evaluation
+================================================================================
+
+.. versionadded:: 5.0.0
+
++-----------------+-------------------------------------------------------------------+
+| **Alias**       | horner                                                            |
++-----------------+-------------------------------------------------------------------+
+| **Domain**      | 2D and 3D                                                         |
++-----------------+-------------------------------------------------------------------+
+| **Input type**  | Geodetic and projected coordinates                                |
++-----------------+-------------------------------------------------------------------+
+| **Output type** | Geodetic and projected coordinates                                |
++-----------------+-------------------------------------------------------------------+
+
+The Horner polynomial evaluation scheme is used for transformations between reference
+frames where one or both are inhomogenous or internally distorted. This will typically
+be reference frames created before the introduction of space geodetic techniques such
+as GPS.
+
+Horner polynomials, or Multiple Regression Equations as they are also known as, have
+their strength in being able to create complicated mappings between coordinate
+reference frames while still being lightweight in both computational cost and disk
+space used.
+
+PROJ implements two versions of the Horner evaluation scheme: Real and complex
+polynomial evaluation. Below both are briefly described. For more details consult
+:cite:`Ruffhead2016` and :cite:`EPSGGuidanceNumber7Part2`.
+
+The polynomial evaluation in real number space is defined by the following equations:
+
+.. math::
+    :label: real_poly
+
+    \Delta X = \sum_{i,j} u_{i,j} U^i V^j
+
+    \Delta Y = \sum_{i,j} v_{i,j} U^i V^j
+
+where
+
+.. math::
+    :label: UV
+
+    U = X_{in} - X_{origin}
+
+    V = Y_{in} - Y_{origin}
+
+and :math:`u_{i,j}` and :math:`v_{i,j}` are coefficients that make up
+the polynomial.
+
+The final coordinates are determined as
+
+.. math::
+    :label: xy_out
+
+    X_{out} = X_{in} + \Delta X
+
+    Y_{out} = Y_{in} + \Delta Y
+
+The inverse transform is the same as the above but requires a different set of
+coefficients.
+
+Evaluation of the complex polynomials are defined by the following equations:
+
+.. math::
+    :label: complex_poly
+
+    \Delta X + i \Delta Y = \sum_{j=1}^n (c_{2j-1} + i c_{2j}) (U + i V)^j
+
+Where :math:`n` is the degree of the polynomial. :math:`U` and :math:`V` are
+defined as in :eq:`UV` and the resulting coordinates are again determined
+by :eq:`xy_out`.
+
+Examples
+################################################################################
+
+Mapping between Danish TC32 and ETRS89/UTM zone 32 using polynomials in real
+number space::
+
+    +proj=horner
+    +ellps=intl
+    +range=500000
+    +fwd_origin=877605.269066,6125810.306769
+    +inv_origin=877605.760036,6125811.281773
+    +deg=4
+    +fwd_v=6.1258112678e+06,9.9999971567e-01,1.5372750011e-10,5.9300860915e-15,2.2609497633e-19,4.3188227445e-05,2.8225130416e-10,7.8740007114e-16,-1.7453997279e-19,1.6877465415e-10,-1.1234649773e-14,-1.7042333358e-18,-7.9303467953e-15,-5.2906832535e-19,3.9984284847e-19
+    +fwd_u=8.7760574982e+05,9.9999752475e-01,2.8817299305e-10,5.5641310680e-15,-1.5544700949e-18,-4.1357045890e-05,4.2106213519e-11,2.8525551629e-14,-1.9107771273e-18,3.3615590093e-10,2.4380247154e-14,-2.0241230315e-18,1.2429019719e-15,5.3886155968e-19,-1.0167505000e-18
+    +inv_v=6.1258103208e+06,1.0000002826e+00,-1.5372762184e-10,-5.9304261011e-15,-2.2612705361e-19,-4.3188331419e-05,-2.8225549995e-10,-7.8529116371e-16,1.7476576773e-19,-1.6875687989e-10,1.1236475299e-14,1.7042518057e-18,7.9300735257e-15,5.2881862699e-19,-3.9990736798e-19
+    +inv_u=8.7760527928e+05,1.0000024735e+00,-2.8817540032e-10,-5.5627059451e-15,1.5543637570e-18,4.1357152105e-05,-4.2114813612e-11,-2.8523713454e-14,1.9109017837e-18,-3.3616407783e-10,-2.4382678126e-14,2.0245020199e-18,-1.2441377565e-15,-5.3885232238e-19,1.0167203661e-18
+
+
+Mapping between Danish System Storebælt and ETRS89/UTM zone 32 using complex
+polynomials::
+
+    +proj=horner
+    +ellps=intl
+    +range=500000
+    +fwd_origin=4.94690026817276e+05,6.13342113183056e+06
+    +inv_origin=6.19480258923588e+05,6.13258568148837e+06
+    +deg=3
+    +fwd_c=6.13258562111350e+06,6.19480105709997e+05,9.99378966275206e-01,-2.82153291753490e-02,-2.27089979140026e-10,-1.77019590701470e-09,1.08522286274070e-14,2.11430298751604e-15
+    +inv_c=6.13342118787027e+06,4.94690181709311e+05,9.99824464710368e-01,2.82279070814774e-02,7.66123542220864e-11,1.78425334628927e-09,-1.05584823306400e-14,-3.32554258683744e-15
+
+Parameters
+################################################################################
+
+Setting up Horner polynomials requires many coefficients being explicitly
+written, even for polynomials of low degree. For this reason it is recommended
+to store the polynomial definitions in an :ref:`init file <init_files>` for
+easier writing and reuse.
+
+Required
+-------------------------------------------------------------------------------
+
+Below is a list of required parameters that can be set for the Horner polynomial
+transformation. As stated above, the transformation takes to forms, either using
+real or complex polynomials. These are divided into separate sections below.
+Parameters from the two sections are mutually exclusive, that is parameters
+describing real and complex polynomials can't be mixed.
+
+.. include:: ../options/ellps.rst
+
+.. option:: +deg=<value>
+
+    Degree of polynomial
+
+.. option:: +fwd_origin=<northing,easting>
+
+    Coordinate of origin for the forward mapping
+
+.. option:: +inv_origin=<northing,easting>
+
+    Coordinate of origin for the inverse mapping
+
+Real polynomials
+..............................................................................
+
+The following parameters has to be set if the transformation consists of
+polynomials in real space. Each parameter takes a comma-separated list of
+coefficients. The number of coefficients is governed by the degree, :math:`d`,
+of the polynomial:
+
+.. math::
+
+    N = \frac{(d + 1)(d + 2)}{2}
+
+.. option:: +fwd_u=<u_11,u_12,...,u_ij,..,u_mn>
+
+    Coefficients for the forward transformation i.e. latitude to northing
+    as described in :eq:`real_poly`.
+
+.. option:: +fwd_v=<v_11,v_12,...,v_ij,..,v_mn>
+
+    Coefficients for the forward transformation i.e. longitude to easting
+    as described in :eq:`real_poly`.
+
+.. option:: +inv_u=<u_11,u_12,...,u_ij,..,u_mn>
+
+    Coefficients for the inverse transformation i.e. latitude to northing
+    as described in :eq:`real_poly`.
+
+.. option:: +inv_v=<v_11,v_12,...,v_ij,..,v_mn>
+
+    Coefficients for the inverse transformation i.e. longitude to easting
+    as described in :eq:`real_poly`.
+
+
+
+Complex polynomials
+..............................................................................
+
+The following parameters has to be set if the transformation consists of
+polynomials in complex space. Each parameter takes a comma-separated list of
+coefficients. The number of coefficients is governed by the degree, :math:`d`,
+of the polynomial:
+
+.. math::
+
+    N = 2d + 2
+
+.. option:: +fwd_c=<c_1,c_2,...,c_N>
+
+    Coefficients for the complex forward transformation
+    as described in :eq:`complex_poly`.
+
+.. option:: +inv_c=<c_1,c_2,...,c_N>
+
+    Coefficients for the complex inverse transformation
+    as described in :eq:`complex_poly`.
+
+Optional
+-------------------------------------------------------------------------------
+
+.. option:: +range=<value>
+
+    Radius of the region of validity.
+
+.. option:: +uneg
+
+    Express latitude as southing. Only applies for complex polynomials.
+
+.. option:: +vneg
+
+    Express longitude as westing. Only applies for complex polynomials.
+
+
+Further reading
+################################################################################
+
+#. `Wikipedia <https://en.wikipedia.org/wiki/Horner%27s_method>`_
diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst
index 121215f7..45c4a80b 100644
--- a/docs/source/operations/transformations/index.rst
+++ b/docs/source/operations/transformations/index.rst
@@ -14,6 +14,7 @@ systems are based on different datums.
    deformation
    geogoffset
    helmert
+   horner
    molodensky
    hgridshift
    vgridshift
diff --git a/docs/source/references.bib b/docs/source/references.bib
index 0fb077ec..1e6929bb 100644
--- a/docs/source/references.bib
+++ b/docs/source/references.bib
@@ -197,6 +197,18 @@
   Url                      = {https://archive.org/details/DTIC_ADA026294}
 }
 
+@article{Ruffhead2016,
+    author = {A. C. Ruffhead},
+    title = {Introduction to multiple regression equations in datum transformations and their reversibility},
+    journal = {Survey Review},
+    volume = {50},
+    number = {358},
+    pages = {82-90},
+    year  = {2016},
+    publisher = {Taylor & Francis},
+    doi = {10.1080/00396265.2016.1244143},
+}
+
 @Book{Snyder1993,
   Title                    = {Flattening the Earth},
   Author                   = {J. P. Snyder},
@@ -247,6 +259,18 @@
   Edition                  = {15th}
 }
 
+@Article{Tobler2017,
+  Author                   = {W. Tobler},
+  Title                    = {A new companion for Mercator},
+  Journal                  = {Cartography and Geographic Information Science},
+  Volume                   = {45},
+  Number                   = {3},
+  Pages                    = {284-285},
+  Publisher                = {Taylor & Francis},
+  Year                     = {2017},
+  Doi                      = {10.1080/15230406.2017.1308837},
+}
+
 @Article{Verey2017,
   Title                    = {Theoretical analysis and practical consequences of adopting a model {ATPOL} grid as a conical projection defining the conversion of plane coordinates to the {WGS} 84 ellipsoid},
   Author                   = {M. Verey},
diff --git a/src/Makefile.am b/src/Makefile.am
index ae5e0a2e..359f2f28 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -53,9 +53,9 @@ libproj_la_SOURCES = \
 	PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c PJ_gall.c PJ_geoc.c \
 	PJ_labrd.c PJ_lsat.c PJ_misrsom.c PJ_merc.c \
 	PJ_mill.c PJ_ocea.c PJ_omerc.c PJ_somerc.c \
-	PJ_tcc.c PJ_tcea.c PJ_times.c PJ_tmerc.c \
+	PJ_tcc.c PJ_tcea.c PJ_times.c PJ_tmerc.c PJ_tobmerc.c \
 	PJ_airy.c PJ_aitoff.c PJ_august.c PJ_bacon.c \
-	PJ_chamb.c PJ_hammer.c PJ_lagrng.c PJ_larr.c \
+	PJ_bertin1953.c PJ_chamb.c PJ_hammer.c PJ_lagrng.c PJ_larr.c \
 	PJ_lask.c PJ_latlong.c PJ_nocol.c PJ_ob_tran.c PJ_oea.c \
 	PJ_tpeqd.c PJ_vandg.c PJ_vandg2.c PJ_vandg4.c \
 	PJ_wag7.c PJ_lcca.c PJ_geos.c proj_etmerc.c \
diff --git a/src/PJ_bertin1953.c b/src/PJ_bertin1953.c
new file mode 100644
index 00000000..5a027da2
--- /dev/null
+++ b/src/PJ_bertin1953.c
@@ -0,0 +1,94 @@
+/*
+  Created by Jacques Bertin in 1953, this projection was the go-to choice
+  of the French cartographic school when they wished to represent phenomena
+  on a global scale.
+  
+  Formula designed by Philippe Rivière, 2017.
+  https://visionscarto.net/bertin-projection-1953
+
+  Port to PROJ by Philippe Rivière, 21 September 2018
+*/
+
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj_internal.h"
+#include "proj.h"
+#include "projects.h"
+
+PROJ_HEAD(bertin1953, "Bertin 1953")
+    "\n\tMisc Sph no inv.";
+
+struct pj_opaque {
+    double cos_delta_phi, sin_delta_phi, cos_delta_gamma, sin_delta_gamma, deltaLambda;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+
+    double fu = 1.4, k = 12., w = 1.68, d;
+
+    /* Rotate */
+    double cosphi, x, y, z, z0;
+    lp.lam += PJ_TORAD(-16.5);
+    cosphi = cos(lp.phi);
+    x = cos(lp.lam) * cosphi;
+    y = sin(lp.lam) * cosphi;
+    z = sin(lp.phi);
+    z0 = z * Q->cos_delta_phi + x * Q->sin_delta_phi;
+    lp.lam = atan2(y * Q->cos_delta_gamma - z0 * Q->sin_delta_gamma,
+       x * Q->cos_delta_phi - z * Q->sin_delta_phi);
+    z0 = z0 * Q->cos_delta_gamma + y * Q->sin_delta_gamma;
+    lp.phi = asin(z0);
+
+    lp.lam = adjlon(lp.lam);
+
+    /* Adjust pre-projection */
+    if (lp.lam + lp.phi < -fu) {
+      d = (lp.lam - lp.phi + 1.6) * (lp.lam + lp.phi + fu) / 8.;
+      lp.lam += d;
+      lp.phi -= 0.8 * d * sin(lp.phi + M_PI / 2.);
+    }
+
+    /* Project with Hammer (1.68,2) */
+    cosphi = cos(lp.phi);
+    d = sqrt(2./(1. + cosphi * cos(lp.lam / 2.)));
+    xy.x = w * d * cosphi * sin(lp.lam / 2.);
+    xy.y = d * sin(lp.phi);
+
+    /* Adjust post-projection */
+    d = (1. - cos(lp.lam * lp.phi)) / k;
+    if (xy.y < 0.) {
+      xy.x *= 1. + d;
+    }
+    if (xy.y > 0.) {
+      xy.x *= 1. + d / 1.5 * xy.x * xy.x;
+    }
+
+    return xy;
+}
+
+
+PJ *PROJECTION(bertin1953) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = Q;
+
+    P->lam0 = 0;
+    P->phi0 = PJ_TORAD(-42.);
+
+    Q->cos_delta_phi = cos(P->phi0);
+    Q->sin_delta_phi = sin(P->phi0);
+    Q->cos_delta_gamma = 1.;
+    Q->sin_delta_gamma = 0.;
+
+    P->es = 0.;
+    P->fwd = s_forward;
+
+    return P;
+}
diff --git a/src/PJ_lcc.c b/src/PJ_lcc.c
index 34ed99cb..96aa3f2a 100644
--- a/src/PJ_lcc.c
+++ b/src/PJ_lcc.c
@@ -5,9 +5,9 @@
 #include "proj_math.h"
 
 PROJ_HEAD(lcc, "Lambert Conformal Conic")
-    "\n\tConic, Sph&Ell\n\tlat_1= and lat_2= or lat_0";
+    "\n\tConic, Sph&Ell\n\tlat_1= and lat_2= or lat_0, k_0=";
 
-# define EPS10  1.e-10
+#define EPS10 1.e-10
 
 struct pj_opaque {
     double phi1;
@@ -15,12 +15,11 @@ struct pj_opaque {
     double n;
     double rho0;
     double c;
-    int    ellips;
 };
 
 
 static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
-    XY xy = {0.0,0.0};
+    XY xy = {0., 0.};
     struct pj_opaque *Q = P->opaque;
     double rho;
 
@@ -31,18 +30,19 @@ static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
         }
         rho = 0.;
     } else {
-        rho = Q->c * (Q->ellips ? pow(pj_tsfn(lp.phi, sin(lp.phi),
-            P->e), Q->n) : pow(tan(M_FORTPI + .5 * lp.phi), -Q->n));
+        rho = Q->c * (P->es != 0. ?
+                      pow(pj_tsfn(lp.phi, sin(lp.phi), P->e), Q->n) :
+                      pow(tan(M_FORTPI + .5 * lp.phi), -Q->n));
     }
     lp.lam *= Q->n;
-    xy.x = P->k0 * (rho * sin( lp.lam) );
-    xy.y = P->k0 * (Q->rho0 - rho * cos(lp.lam) );
+    xy.x = P->k0 * (rho * sin(lp.lam));
+    xy.y = P->k0 * (Q->rho0 - rho * cos(lp.lam));
     return xy;
 }
 
 
 static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
-    LP lp = {0.0,0.0};
+    LP lp = {0., 0.};
     struct pj_opaque *Q = P->opaque;
     double rho;
 
@@ -51,13 +51,13 @@ static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
 
     xy.y = Q->rho0 - xy.y;
     rho = hypot(xy.x, xy.y);
-    if (rho != 0.0) {
+    if (rho != 0.) {
         if (Q->n < 0.) {
             rho = -rho;
             xy.x = -xy.x;
             xy.y = -xy.y;
         }
-        if (Q->ellips) {
+        if (P->es != 0.) {
             lp.phi = pj_phi2(P->ctx, pow(rho / Q->c, 1./Q->n), P->e);
             if (lp.phi == HUGE_VAL) {
                 proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
@@ -78,13 +78,12 @@ static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
 PJ *PROJECTION(lcc) {
     double cosphi, sinphi;
     int secant;
-    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    struct pj_opaque *Q = pj_calloc(1, sizeof (struct pj_opaque));
 
-    if (0==Q)
-        return pj_default_destructor (P, ENOMEM);
+    if (0 == Q)
+        return pj_default_destructor(P, ENOMEM);
     P->opaque = Q;
 
-
     Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
     if (pj_param(P->ctx, P->params, "tlat_2").i)
         Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
@@ -99,10 +98,9 @@ PJ *PROJECTION(lcc) {
     Q->n = sinphi = sin(Q->phi1);
     cosphi = cos(Q->phi1);
     secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
-    if( (Q->ellips = (P->es != 0.)) ) {
+    if (P->es != 0.) {
         double ml1, m1;
 
-        P->e = sqrt(P->es);
         m1 = pj_msfn(sinphi, cosphi, P->es);
         ml1 = pj_tsfn(Q->phi1, sinphi, P->e);
         if (secant) { /* secant cone */
@@ -128,4 +126,3 @@ PJ *PROJECTION(lcc) {
 
     return P;
 }
-
diff --git a/src/PJ_tobmerc.c b/src/PJ_tobmerc.c
new file mode 100644
index 00000000..9c939f0b
--- /dev/null
+++ b/src/PJ_tobmerc.c
@@ -0,0 +1,51 @@
+#define PJ_LIB__
+
+#include <float.h>
+#include <math.h>
+
+#include "proj_internal.h"
+#include "proj.h"
+#include "proj_math.h"
+#include "projects.h"
+
+PROJ_HEAD(tobmerc, "Tobler-Mercator") "\n\tCyl, Sph";
+
+#define EPS10 1.e-10
+static double logtanpfpim1(double x) {       /* log(tan(x/2 + M_FORTPI)) */
+    if (fabs(x) <= DBL_EPSILON) {
+        /* tan(M_FORTPI + .5 * x) can be approximated by  1.0 + x */
+        return log1p(x);
+    }
+    return log(tan(M_FORTPI + .5 * x));
+}
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0, 0.0};
+    double cosphi;
+
+    if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
+        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+        return xy;
+    }
+
+    cosphi = cos(lp.phi);
+    xy.x = P->k0 * lp.lam * cosphi * cosphi;
+    xy.y = P->k0 * logtanpfpim1(lp.phi);
+    return xy;
+}
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0, 0.0};
+    double cosphi;
+
+    lp.phi = atan(sinh(xy.y / P->k0));
+    cosphi = cos(lp.phi);
+    lp.lam = xy.x / P->k0 / (cosphi * cosphi);
+    return lp;
+}
+
+PJ *PROJECTION(tobmerc) {
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+    return P;
+}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index a359369e..6a287a43 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -46,6 +46,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_august.c
         PJ_axisswap.c
         PJ_bacon.c
+        PJ_bertin1953.c
         PJ_bipc.c
         PJ_boggs.c
         PJ_bonne.c
@@ -140,6 +141,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_tcea.c
         PJ_times.c
         PJ_tmerc.c
+        PJ_tobmerc.c
         PJ_tpeqd.c
         PJ_unitconvert.c
         PJ_urm5.c
diff --git a/src/pj_gridinfo.c b/src/pj_gridinfo.c
index f201f39e..de0e8d31 100644
--- a/src/pj_gridinfo.c
+++ b/src/pj_gridinfo.c
@@ -660,7 +660,7 @@ static int pj_gridinfo_init_ntv2( projCtx ctx, PAFile fid, PJ_GRIDINFO *gilist )
 static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi )
 
 {
-    unsigned char header[176];
+    unsigned char header[192]; /* 12 records of 16 bytes */
     struct CTABLE *ct;
     LP		ur;
 
@@ -731,7 +731,7 @@ static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi )
     ct->cvs = NULL;
 
     gi->ct = ct;
-    gi->grid_offset = pj_ctx_ftell( ctx, fid );
+    gi->grid_offset = (long) sizeof(header);
     gi->format = "ntv1";
 
     return 1;
diff --git a/src/pj_list.h b/src/pj_list.h
index 72a9af7f..8c72dd1f 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -15,6 +15,7 @@ PROJ_HEAD(apian, "Apian Globular I")
 PROJ_HEAD(august, "August Epicycloidal")
 PROJ_HEAD(axisswap, "Axis ordering")
 PROJ_HEAD(bacon, "Bacon Globular")
+PROJ_HEAD(bertin1953, "Bertin 1953")
 PROJ_HEAD(bipc, "Bipolar conic of western hemisphere")
 PROJ_HEAD(boggs, "Boggs Eumorphic")
 PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)")
@@ -139,6 +140,7 @@ PROJ_HEAD(tcea, "Transverse Cylindrical Equal Area")
 PROJ_HEAD(times, "Times Projection")
 PROJ_HEAD(tissot, "Tissot Conic")
 PROJ_HEAD(tmerc, "Transverse Mercator")
+PROJ_HEAD(tobmerc, "Tobler-Mercator")
 PROJ_HEAD(tpeqd, "Two Point Equidistant")
 PROJ_HEAD(tpers, "Tilted perspective")
 PROJ_HEAD(unitconvert, "Unit conversion")
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index 1ee968f4..f4e15a77 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -2505,11 +2505,10 @@ expect  -217928.275907355 112144.329220142
 accept  -2 -1
 expect  -217928.275907355 -112144.329220142
 
-
 ===============================================================================
 Lambert Conformal Conic
 	Conic, Sph&Ell
-	lat_1= and lat_2= or lat_0
+	lat_1= and lat_2= or lat_0, k_0=
 ===============================================================================
 
 -------------------------------------------------------------------------------
@@ -2535,6 +2534,94 @@ expect  -0.001796359 0.000904232
 accept  -200 -100
 expect  -0.001796358 -0.000904233
 
+-------------------------------------------------------------------------------
+Tests with +k_0 (Lambert Conformal Conic 2SP Michigan).
+-------------------------------------------------------------------------------
+operation +proj=lcc +ellps=GRS80 +lat_1=0.5 +lat_2=2 +k_0=1.0000382
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  2 1
+expect  222596.942614366 110664.761103214
+accept  2 -1
+expect  222765.389013083 -110537.020013748
+accept  -2 1
+expect  -222596.942614366 110664.761103214
+accept  -2 -1
+expect  -222765.389013083 -110537.020013748
+
+direction inverse
+accept  200 100
+expect  0.001796291 0.000904198
+accept  200 -100
+expect  0.001796290 -0.000904199
+accept  -200 100
+expect  -0.001796291 0.000904198
+accept  -200 -100
+expect  -0.001796290 -0.000904199
+
+
+-------------------------------------------------------------------------------
+Test various corner cases, spherical projection and one standard
+parallel to improve test coverage.
+-------------------------------------------------------------------------------
+operation +proj=lcc +ellps=GRS80 +lat_1=0.5 +lat_2=2
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept  0 90
+expect  0 292411117.537843227
+accept 0 0
+expect 0 0
+
+-------------------------------------------------------------------------------
+operation +proj=lcc +ellps=sphere +lat_1=30 +lat_2=40
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept 1 2
+expect 129391.909521100 262101.674176860
+accept 0 0
+expect 0 0
+
+direction inverse
+
+accept 129391.909521100 262101.674176860
+expect 1 2
+accept 0 0
+expect 0 0
+
+-------------------------------------------------------------------------------
+operation +proj=lcc +ellps=GRS80 +lat_1=30
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept 1 2
+expect 131833.493971117 265456.213515346
+accept 1 -2
+expect 137536.205750651 -269686.591917190
+accept -1 2
+expect -131833.493971117 265456.213515346
+accept -1 -2
+expect -137536.205750651 -269686.591917190
+
+-------------------------------------------------------------------------------
+operation +proj=lcc +ellps=sphere +lat_1=30
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept 1 2
+expect 131824.206082557 267239.875053699
+accept -1 -2
+expect -137565.475967350 -271546.945608449
+accept 1 -2
+expect 137565.475967350 -271546.945608449
+accept -1 2
+expect -131824.206082557 267239.875053699
+
+direction inverse
+
+accept 131824.206082557	267239.875053699
+expect 1 2
 
 ===============================================================================
 Lambert Conformal Conic Alternative
@@ -4810,6 +4897,72 @@ expect  -0.001790493 0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+===============================================================================
+Tobler-Mercator
+	Cyl, Sph
+===============================================================================
+
+-------------------------------------------------------------------------------
+operation +proj=tobmerc +ellps=sphere
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  2 1
+expect  222322.011656333081 111200.520030584055
+accept  2 -1
+expect  222322.011656333081 -111200.520030584463
+accept  -2 1
+expect  -222322.011656333081 111200.520030584055
+accept  -2 -1
+expect  -222322.011656333081 -111200.520030584463
+
+accept  2 75
+expect  14897.288383530242 12917766.123520378023
+accept  2 80
+expect  6705.871450149517 15521316.299485698342
+accept  -2 -75
+expect  -14897.288383530242 -12917766.123520381749
+accept  -2 -80
+expect  -6705.871450149517 -15521316.299485692754
+
+direction inverse
+accept  200 100
+expect  0.001798644059 0.000899322029
+accept  200 -100
+expect  0.001798644059 -0.000899322029
+accept  -200 100
+expect  -0.001798644059 0.000899322029
+accept  -200 -100
+expect  -0.001798644059 -0.000899322029
+
+accept  14897.288383530242 12917766.123520378023
+expect  2 75
+accept  6705.871450149517 15521316.299485698342
+expect  2 80
+accept  -14897.288383530242 -12917766.123520381749
+expect  -2 -75
+accept  -6705.871450149517 -15521316.299485692754
+expect  -2 -80
+
+-------------------------------------------------------------------------------
+operation +proj=tobmerc +R=1
+-------------------------------------------------------------------------------
+Test the numerical stability of the inverse spherical Tobler-Mercator
+-------------------------------------------------------------------------------
+tolerance   1e-15 m
+accept  0   1e-15
+expect  0   1e-15
+
+-------------------------------------------------------------------------------
+operation +proj=tobmerc +ellps=sphere
+-------------------------------------------------------------------------------
+Test expected failure at the poles:
+-------------------------------------------------------------------------------
+accept  0 90
+expect  failure errno tolerance_condition
+
+accept  0 -90
+expect  failure errno tolerance_condition
+
 
 ===============================================================================
 Two Point Equidistant
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 10feea70..6ac0c70c 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -59,6 +59,41 @@ roundtrip 100    1 m
 -------------------------------------------------------------------------------
 
 
+-------------------------------------------------------------------------------
+Tests for PJ_bertin1953.c
+-------------------------------------------------------------------------------
+operation proj=bertin1953 +R=1
+-------------------------------------------------------------------------------
+accept    0 0
+expect    -0.260206554508 -0.685226058142
+
+accept    16.5 42
+expect    0.0 0.0
+
+accept    -180 90
+expect    0.0 0.813473286152
+
+accept    0 90
+expect    0.0 0.813473286152
+
+accept    10 -35
+expect    -0.138495501548 -1.221408328101
+
+accept    -70 -35
+expect    -1.504967424950 -0.522846035499
+
+accept    80 7
+expect    0.929377425352 -0.215443296201
+
+accept    128 35
+expect    0.937431733950 0.700084268525
+
+accept    170 -41
+expect    2.162845830414 -0.046534568425
+
+-------------------------------------------------------------------------------
+
+
 
 
 
@@ -243,6 +278,23 @@ expect    failure errno no_args
 
 
 
+-------------------------------------------------------------------------------
+Tests for LCC 2SP Michigan (from PJ_lcc.c)
+-------------------------------------------------------------------------------
+This test is taken from EPSG guidance note 7-2 (version 54, August 2018,
+page 25)
+-------------------------------------------------------------------------------
+operation +proj=lcc +ellps=clrk66 +lat_1=44d11'N +lat_2=45d42'N +x_0=609601.2192 +lon_0=84d20'W +lat_0=43d19'N +k_0=1.0000382 +units=us-ft
+-------------------------------------------------------------------------------
+tolerance 5 mm
+accept 83d10'W 43d45'N
+expect 2308335.75 160210.48
+
+direction inverse
+accept 2308335.75 160210.48
+expect 83d10'W 43d45'N
+-------------------------------------------------------------------------------
+
 
 -------------------------------------------------------------------------------
 A number of tests from PJ_helmert.c
diff --git a/test/old/td_out.dist b/test/old/td_out.dist
index 060d14d2..f6b2a219 100644
--- a/test/old/td_out.dist
+++ b/test/old/td_out.dist
@@ -1,6 +1,6 @@
 ##############################################################
 1st through ntv1, 2nd through conus
-111d00'00.000"W 44d00'00.000"N 0.0	111d0'3.085"W	43d59'59.756"N 0.000
+111d00'00.000"W 44d00'00.000"N 0.0	111d0'3.208"W	43d59'59.732"N 0.000
 111d00'00.000"W 39d00'00.000"N 0.0	111d0'2.604"W	38d59'59.912"N 0.000
 ##############################################################
 As above, but without ntv1 everything goes through conus file.