Merge remote-tracking branch 'osgeo/master' into remove-t_obs

63 files changed, 1790 insertions, 1932 deletions

diff --git a/cmake/ProjUtilities.cmake b/cmake/ProjUtilities.cmake
index 6969a664..72373719 100644
--- a/cmake/ProjUtilities.cmake
+++ b/cmake/ProjUtilities.cmake
@@ -111,57 +111,4 @@ function(proj_target_output_name TARGET_NAME OUTPUT_NAME)
   endif()
 
   set(${OUTPUT_NAME} ${TARGET_NAME}${SUFFIX} PARENT_SCOPE)
-endfunction()
-
-
-#
-# conversion from lla name to lla convert name ( without lla extension)
-#
-
-function(proj_lla_output_name LLA_INPUT_NAME LLA_OUTPUT_NAME  )
-    get_filename_component(filename ${LLA_INPUT_NAME} NAME_WE)
-    get_filename_component(pathname ${LLA_INPUT_NAME} PATH)
-    set(${LLA_OUTPUT_NAME} ${pathname}/${filename} PARENT_SCOPE)
-        set(${LLA_OUTPUT_NAME} ${pathname}/${filename} PARENT_SCOPE)
-endfunction()
-
-function(proj_lla_target_name LLA_INPUT_NAME  LLA_TARGET )
-    get_filename_component(filename ${LLA_INPUT_NAME} NAME_WE)
-    set(${LLA_TARGET} ${filename} PARENT_SCOPE)
-endfunction()
-
-#
-# in place conversion of lla file to gsb 
-#
-
-function(proj_convert_grid_lla2gsb GRID_DIRECTORY) 
-    set(NAD2BIN_DIR ${CMAKE_RUNTIME_OUTPUT_DIRECTORY})
-    set(NAD2BIN_PATH ${NAD2BIN_DIR}/nad2bin${CMAKE_EXECUTABLE_SUFFIX})
-    file(TO_NATIVE_PATH ${NAD2BIN_PATH} NAD2BIN_EXE)
-    file(GLOB LLA_FILES  ${${GRID_DIRECTORY}}/*.lla)
-    foreach(LLA ${LLA_FILES} )
-        proj_lla_output_name(${LLA} DEST_FILE)
-        file(TO_NATIVE_PATH ${DEST_FILE} DEST)
-        proj_lla_target_name(${LLA} LLA_TARGET)
-        if(NOT EXISTS ${DEST})
-            add_custom_target( ${LLA_TARGET} ALL
-               COMMAND ${NAD2BIN_EXE} ${DEST} "<" ${LLA}
-               DEPENDS nad2bin )
-        endif(NOT EXISTS ${DEST})
-    endforeach(LLA)
-endfunction()
-
-#
-# add lla output list to an existing file list
-#
-
-function(proj_append_lla_output_file LLA_INPUT_FILE  FILE_LIST)
-     set(LIST_OUT ${${FILE_LIST}})
-     foreach(LLA ${${LLA_INPUT_FILE}} )
-        proj_lla_output_name(${LLA} DEST_FILE)
-        file(TO_NATIVE_PATH ${DEST_FILE} DEST)
-        set(LIST_OUT ${LIST_OUT} ${DEST_FILE} )
-    endforeach(LLA ${LLA_INPUT_FILE})
-    set(${FILE_LIST} ${LIST_OUT} PARENT_SCOPE)
-endfunction()
-
+endfunction()
\ No newline at end of file
diff --git a/configure.ac b/configure.ac
index 95d587c6..1135b4a4 100644
--- a/configure.ac
+++ b/configure.ac
@@ -309,37 +309,9 @@ dnl ---------------------------------------------------------------------------
 dnl Check for sqlite3 library and binary
 dnl ---------------------------------------------------------------------------
 
-AC_ARG_WITH(sqlite3-include,
-    [  --with-sqlite3-include=ARG sqlite3 include directory],,,)
-
-if test "x$with_sqlite3_include" = "x" ; then
-SQLITE3_FLAGS=""
-if test -z "`uname | grep Darwin`" ; then
-    AC_CHECK_HEADERS([sqlite3.h])
-    if test "$ac_cv_header_sqlite3_h" = "no" ; then
-        AC_MSG_ERROR([Please install sqlite3 development package.])
-    fi
-fi
-else
-    SQLITE3_FLAGS="-I$with_sqlite3_include"
-fi
-AC_SUBST(SQLITE3_FLAGS,$SQLITE3_FLAGS)
-
-AC_ARG_WITH(sqlite3-ldflags,
-      [  --with-sqlite3-ldflags=ARG   sqlite3 link flags (i.e. -L/xxx -lsqlite3)],,,)
-
-SQLITE3_LDFLAGS=""
-if test "x$with_sqlite3_ldflags" = "x" ; then
-    AC_CHECK_LIB(sqlite3,sqlite3_open,SQLITE3_LIB_FOUND=yes,SQLITE3_LIB_FOUND=no,)
-    if test x"$SQLITE3_LIB_FOUND" != x"yes"; then
-        AC_MSG_ERROR([Please install sqlite3 development package.])
-    fi
-    SQLITE3_LDFLAGS="-lsqlite3"
-else
-    SQLITE3_LDFLAGS="$with_sqlite3_ldflags"
-fi
-AC_SUBST(SQLITE3_LDFLAGS,$SQLITE3_LDFLAGS)
-
+PKG_CHECK_MODULES([SQLITE3], [sqlite3 >= 3.7])
+AC_SUBST(SQLITE3_CFLAGS,$SQLITE3_CFLAGS)
+AC_SUBST(SQLITE3_LIBS,$SQLITE3_LIBS)
 
 AC_CHECK_PROG(SQLITE3_CHECK,sqlite3,yes)
 if test x"$SQLITE3_CHECK" != x"yes" ; then
@@ -347,19 +319,40 @@ if test x"$SQLITE3_CHECK" != x"yes" ; then
 fi
 
 dnl ---------------------------------------------------------------------------
+dnl Check for external Google Test
+dnl ---------------------------------------------------------------------------
+
+AC_ARG_WITH(external-gtest,
+            AS_HELP_STRING([--with-external-gtest],
+               [Whether to use external Google Test]),,)
+
+if test "x$with_external_gtest" == "xyes" ; then
+    AC_MSG_RESULT([using external GTest.])
+    PKG_CHECK_MODULES([GTEST], [gtest >= 1.8.0])
+else
+    AC_MSG_RESULT([using internal GTest.])
+    GTEST_CFLAGS="-I\$(top_srcdir)/test/googletest/include"
+    GTEST_LIBS="\$(top_builddir)/test/googletest/libgtest.la"
+fi
+AM_CONDITIONAL(USE_EXTERNAL_GTEST, [test "x$with_external_gtest" == "xyes"])
+AC_SUBST(GTEST_CFLAGS,$GTEST_CFLAGS)
+AC_SUBST(GTEST_LIBS,$GTEST_LIBS)
+
+dnl ---------------------------------------------------------------------------
 dnl Generate files
 dnl ---------------------------------------------------------------------------
 
 AC_CONFIG_FILES([Makefile cmake/Makefile src/Makefile include/Makefile include/proj/Makefile include/proj/internal/Makefile
                  test/Makefile test/cli/Makefile test/gie/Makefile test/gigs/Makefile test/unit/Makefile
-                 test/googletest/Makefile test/googletest/include/Makefile
-                 test/googletest/include/gtest/Makefile
-                 test/googletest/include/gtest/internal/Makefile
-                 test/googletest/include/gtest/internal/custom/Makefile
-                 test/googletest/src/Makefile
                  man/Makefile man/man1/Makefile man/man3/Makefile data/Makefile
                  jniwrap/Makefile jniwrap/org.osgeo.proj/Makefile jniwrap/org.osgeo.proj/org/Makefile jniwrap/org.osgeo.proj/org/proj4/Makefile])
-AC_CONFIG_FILES([data/install], [chmod +x data/install])
+if ! test "x$with_external_gtest" = "xyes" ; then
+    AC_CONFIG_FILES([test/googletest/Makefile test/googletest/include/Makefile
+                     test/googletest/include/gtest/Makefile
+                     test/googletest/include/gtest/internal/Makefile
+                     test/googletest/include/gtest/internal/custom/Makefile
+                     test/googletest/src/Makefile])
+fi
 
 AC_CONFIG_FILES([proj.pc])
 
diff --git a/data/CMakeLists.txt b/data/CMakeLists.txt
index 859f9a60..8d3a06b2 100644
--- a/data/CMakeLists.txt
+++ b/data/CMakeLists.txt
@@ -22,17 +22,6 @@ file(GLOB GTX_FILES  *.gtx)
 set(GRIDSHIFT_FILES ${GSB_FILES} ${GTX_FILES})
 set(GRIDSHIFT_FILES ${GRIDSHIFT_FILES}
                      )
-option(CONVERT_DATA "convert some ascii file to binary file for use in proj4" OFF)
-if(CONVERT_DATA)
-    if (BUILD_NAD2BIN)
-       set(LLA_GRID_DIR ${CMAKE_CURRENT_SOURCE_DIR})
-       proj_convert_grid_lla2gsb(LLA_GRID_DIR)
-       file(GLOB LLA_FILES  *.lla)
-       proj_append_lla_output_file(LLA_FILES GRIDSHIFT_FILES)
-    else()
-       message(ERROR " you need to compile nad2bin exe in order to convert data file" )
-    endif(BUILD_NAD2BIN)
-endif(CONVERT_DATA)
 
 set(ALL_SQL_IN "${CMAKE_CURRENT_SOURCE_DIR}/all.sql.in")
 set(PROJ_DB "${CMAKE_CURRENT_SOURCE_DIR}/proj.db")
diff --git a/data/Makefile.am b/data/Makefile.am
index 3e98e59f..8234c034 100644
--- a/data/Makefile.am
+++ b/data/Makefile.am
@@ -1,7 +1,7 @@
 DATAPATH = $(top_srcdir)/data
 
 pkgdata_DATA = GL27 nad.lst nad27 nad83 world other.extra \
-		CH \
+		CH null \
 		ITRF2000 ITRF2008 ITRF2014 proj.db
 
 SQL_ORDERED_LIST = sql/begin.sql \
@@ -40,34 +40,12 @@ EXTRA_DIST = GL27 nad.lst nad27 nad83 \
 		world other.extra \
 		CH \
 		ITRF2000 ITRF2008 ITRF2014 \
-		CMakeLists.txt tests/test_nodata.gtx null.lla \
+		CMakeLists.txt tests/test_nodata.gtx null \
 		generate_all_sql_in.cmake sql_filelist.cmake \
 		$(SQL_ORDERED_LIST)
 
-process-nad2bin:
-	@if [ -f $(DATAPATH)/null.lla -a ! -f null ] || [ -f $(DATAPATH)/conus.lla -a ! -f conus ] ; then \
-	  for x in $(DATAPATH)/*.lla ; do \
-	    DEST=`basename $$x .lla`; \
-	    echo "../src/nad2bin$(EXEEXT) $$DEST < $$x"; \
-	    ../src/nad2bin$(EXEEXT) $$DEST < $$x ; \
-	  done; \
-	fi
-
-install-data-local: process-nad2bin
+install-data-local:
 	$(mkinstalldirs) $(DESTDIR)$(pkgdatadir)
-	@if [ -f $(DATAPATH)/null.lla -a ! -f $(DATAPATH)/null ] || [ -f $(DATAPATH)/conus.lla -a ! -f $(DATAPATH)/conus ] ; then \
-	  for x in $(DATAPATH)/*.lla ; do \
-	    DEST=$(DESTDIR)$(pkgdatadir)/`basename $$x .lla`; \
-	    echo "../src/nad2bin$(EXEEXT) $$DEST < $$x"; \
-	    ../src/nad2bin$(EXEEXT) $$DEST < $$x ; \
-	    if [ $$? -ne 0 ] ; then \
-		echo "nad2bin$(EXEEXT) failed to install $$DEST"; \
-		exit 1; \
-	    fi; \
-	  done; \
-	else \
-	  echo "nad2nad NADCON source files not present"; \
-	fi
 	@for gridfile in $(DATAPATH)/*.gsb $(DATAPATH)/*.gtx $(DATAPATH)/ntv1_can.dat dummy \
 	                 $(DATAPATH)/alaska $(DATAPATH)/conus $(DATAPATH)/hawaii $(DATAPATH)/null \
 	                 $(DATAPATH)/prvi $(DATAPATH)/stgeorge $(DATAPATH)/stlrnc $(DATAPATH)/stpaul \
@@ -97,9 +75,9 @@ proj.db: $(DATAPATH)/sql/*.sql
 		exit 1; \
 	 fi
 
-	
+
 # For out-of-tree builds, link all file of the source data dir to the generated data
-check-local: process-nad2bin
+check-local:
 	@if [ ! -f GL27 ]; then \
 		for x in $(DATAPATH)/*; do \
 			ln -sf $$x .; \
diff --git a/data/README b/data/README
index 1b18dd2f..ae690729 100644
--- a/data/README
+++ b/data/README
@@ -1,5 +1,4 @@
-The files in this directory are support data for PROJ programs
-as well as test scripts for testing proj installation.
+The files in this directory are support data for PROJ programs.
 
 File Contents:
 
@@ -28,7 +27,7 @@ nad.lst --- Reference list of SPCS States and NGS datum identifiers
 
 IGNF    --- Translation of French Mapping Agency codes into PROJ via init= mechanism.
 
-install     --- Unix shell (sh) script to install nad2nad file and proj
-				init files.  Read beginning for usage instructions.
+install     --- Unix shell (sh) script to install proj init files.
+				Read beginning for usage instructions.
 
 Additional data files are available in https://github.com/OSGeo/proj-datumgrid/
diff --git a/data/install.in b/data/install.in
deleted file mode 100755
index d60dc725..00000000
--- a/data/install.in
+++ /dev/null
@@ -1,53 +0,0 @@
-:
-# SCCSID @(#)install.in	4.8 95/09/24 GIE REL
-#
-# Installation script all PROJ.4 system reference files.
-#
-# This script must be executed after compilation of proj library.
-#
-# Execute as:
-#
-#	install
-#
-# ------------------------------------------------------------
-prefix=/usr/local
-#
-lib=${prefix}/lib/proj.4
-# Copy in "init" files
-for x in proj_def.dat nad27 nad83 world GL27
-do
-	cp $x ${lib}/$x
-	if [ $? -ne 0 ] ; then
-		echo "init copying failed for file $x"
-	else
-		echo "file $x installed"
-	fi
-done
-# just test for conus and assume others are there
-if [ -f conus.lla ] ; then
-# install nad2nad NADCON files
-	echo "installing ${lib}/nad2783"
-	lib=${lib}/nad2783
-	if [ ! -d $lib ] ; then
-		mkdir $lib
-		if [ $? -ne 0 ] ; then
-			echo "failed to make dir $lib"
-			echo "nad2nad lib install stopped"
-			exit 1
-		fi
-	fi
-# Process ascii grid into nad2nad compatible files.
-	for x in *.lla ; do
-		dest=`expr $x : '\(.*\)\.lla'`
-		cat < $x | ../src/nad2bin ${lib}/${dest}
-		if [ $? -ne 0 ] ; then
-			echo "nad2nad file install failed on file ${lib}/${dest}"
-			exit 1
-		fi
-		echo "${lib}/${dest} installed"
-	done
-else
-	echo "nad2nad NADCON files not present"
-fi
-echo "normal completion"
-exit 0
diff --git a/data/null b/data/null
new file mode 100644
index 00000000..8418408b
--- /dev/null
+++ b/data/null
diff --git a/data/null.lla b/data/null.lla
deleted file mode 100644
index eb25285a..00000000
--- a/data/null.lla
+++ /dev/null
@@ -1,5 +0,0 @@
-Null transform for whole world
-   3   3   1  -180.00000   180.00000   -90.00000   90.00000
-0: 0 0 0 0 0 0 
-1: 0 0 0 0 0 0 
-2: 0 0 0 0 0 0 
diff --git a/docs/source/operations/projections/etmerc.rst b/docs/source/operations/projections/etmerc.rst
deleted file mode 100644
index 2a007130..00000000
--- a/docs/source/operations/projections/etmerc.rst
+++ /dev/null
@@ -1,32 +0,0 @@
-.. _etmerc:
-
-********************************************************************************
-Extended Transverse Mercator
-********************************************************************************
-
-.. figure:: ./images/etmerc.png
-   :width: 500 px
-   :align: center
-   :alt:   Extended Transverse Mercator
-
-   proj-string: ``+proj=etmerc +lon_0=-40``
-
-Parameters
-################################################################################
-
-.. note:: All parameters for the projection are optional.
-
-.. include:: ../options/lon_0.rst
-
-.. include:: ../options/lat_0.rst
-
-.. include:: ../options/ellps.rst
-
-.. include:: ../options/R.rst
-
-.. include:: ../options/k_0.rst
-
-.. include:: ../options/x_0.rst
-
-.. include:: ../options/y_0.rst
-
diff --git a/docs/source/operations/projections/images/etmerc.png b/docs/source/operations/projections/images/etmerc.png
deleted file mode 100644
index 35bce083..00000000
--- a/docs/source/operations/projections/images/etmerc.png
+++ /dev/null
diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst
index 572deb1a..e6ec9076 100644
--- a/docs/source/operations/projections/index.rst
+++ b/docs/source/operations/projections/index.rst
@@ -43,7 +43,6 @@ Projections map the spherical 3D space to a flat 2D space.
    eqdc
    eqearth
    euler
-   etmerc
    fahey
    fouc
    fouc_s
diff --git a/docs/source/operations/projections/tmerc.rst b/docs/source/operations/projections/tmerc.rst
index 6d1c145d..fd600001 100644
--- a/docs/source/operations/projections/tmerc.rst
+++ b/docs/source/operations/projections/tmerc.rst
@@ -68,13 +68,19 @@ Example using Gauss-Kruger on Germany area (aka EPSG:31467) ::
 Example using Gauss Boaga on Italy area (EPSG:3004) ::
 
     $ echo 15 42 | proj +proj=tmerc +lat_0=0 +lon_0=15 +k_0=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m +no_defs
-    2520000.00	4649858.60 
+    2520000.00	4649858.60
 
 Parameters
 ################################################################################
 
 .. note:: All parameters for the projection are optional.
 
+.. option:: +approx
+
+    .. versionadded:: 6.0.0
+
+    Use faster, less accurate algorithm for the Transverse Mercator.
+
 .. include:: ../options/lon_0.rst
 
 .. include:: ../options/lat_0.rst
diff --git a/docs/source/operations/projections/utm.rst b/docs/source/operations/projections/utm.rst
index ac65aa1a..8759f4eb 100644
--- a/docs/source/operations/projections/utm.rst
+++ b/docs/source/operations/projections/utm.rst
@@ -66,6 +66,12 @@ Required
 
     Add this flag when using the UTM on the southern hemisphere.
 
+.. option:: +approx
+
+    .. versionadded:: 6.0.0
+
+    Use faster, less accurate algorithm for the Transverse Mercator.
+
 Optional
 -------------------------------------------------------------------------------
 
diff --git a/docs/source/usage/transformation.rst b/docs/source/usage/transformation.rst
index bc5c3b1b..da856a01 100644
--- a/docs/source/usage/transformation.rst
+++ b/docs/source/usage/transformation.rst
@@ -210,8 +210,7 @@ at each grid location. Actually grid shifts are normally computed based on an
 interpolation between the containing four grid points.
 
 PROJ supports use of grid files for shifting between various reference frames.
-The grid shift table formats are ctable (the binary format produced by the PROJ
-``nad2bin`` program), NTv1 (the old Canadian format), and NTv2 (``.gsb`` - the new
+The grid shift table formats are ctable, NTv1 (the old Canadian format), and NTv2 (``.gsb`` - the new
 Canadian and Australian format).
 
 The text in this section is based on the *cs2cs* framework. Gridshifting is off
diff --git a/include/proj/common.hpp b/include/proj/common.hpp
index ad8893f7..d0ecf5d6 100644
--- a/include/proj/common.hpp
+++ b/include/proj/common.hpp
@@ -366,7 +366,8 @@ using ObjectDomainNNPtr = util::nn<ObjectDomainPtr>;
  *
  * \remark Implements ObjectDomain from \ref ISO_19111_2018
  */
-class ObjectDomain : public util::BaseObject, public util::IComparable {
+class PROJ_GCC_DLL ObjectDomain : public util::BaseObject,
+                                  public util::IComparable {
   public:
     //! @cond Doxygen_Suppress
     PROJ_DLL ~ObjectDomain() override;
@@ -420,7 +421,7 @@ using ObjectUsageNNPtr = util::nn<ObjectUsagePtr>;
  *
  * \remark Implements ObjectUsage from \ref ISO_19111_2018
  */
-class ObjectUsage : public IdentifiedObject {
+class PROJ_GCC_DLL ObjectUsage : public IdentifiedObject {
   public:
     //! @cond Doxygen_Suppress
     PROJ_DLL ~ObjectUsage() override;
diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index 4120d707..091efcb7 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -365,7 +365,7 @@ class PROJ_GCC_DLL PROJStringFormatter {
     PROJ_DLL ~PROJStringFormatter();
     //! @endcond
 
-    PROJ_DLL void setUseETMercForTMerc(bool flag);
+    PROJ_DLL void setUseApproxTMerc(bool flag);
 
     PROJ_DLL const std::string &toString() const;
 
@@ -378,7 +378,7 @@ class PROJ_GCC_DLL PROJStringFormatter {
     PROJ_DLL void startInversion();
     PROJ_DLL void stopInversion();
     PROJ_INTERNAL bool isInverted() const;
-    PROJ_INTERNAL bool getUseETMercForTMerc(bool &settingSetOut) const;
+    PROJ_INTERNAL bool getUseApproxTMerc() const;
     PROJ_INTERNAL void setCoordinateOperationOptimizations(bool enable);
 
     PROJ_DLL void
@@ -888,6 +888,38 @@ class PROJ_GCC_DLL AuthorityFactory {
 
     PROJ_DLL std::string getDescriptionText(const std::string &code) const;
 
+    /** CRS information */
+    struct CRSInfo {
+        /** Authority name */
+        std::string authName{};
+        /** Code */
+        std::string code{};
+        /** Name */
+        std::string name{};
+        /** Type */
+        ObjectType type{ObjectType::CRS};
+        /** Whether the object is deprecated */
+        bool deprecated{};
+        /** Whereas the west_lon_degree, south_lat_degree, east_lon_degree and
+        * north_lat_degree fields are valid. */
+        bool bbox_valid{};
+        /** Western-most longitude of the area of use, in degrees. */
+        double west_lon_degree{};
+        /** Southern-most latitude of the area of use, in degrees. */
+        double south_lat_degree{};
+        /** Eastern-most longitude of the area of use, in degrees. */
+        double east_lon_degree{};
+        /** Northern-most latitude of the area of use, in degrees. */
+        double north_lat_degree{};
+        /** Name of the area of use. */
+        std::string areaName{};
+        /** Name of the projection method for a projected CRS. Might be empty
+         * even for projected CRS in some cases. */
+        std::string projectionMethodName{};
+    };
+
+    PROJ_DLL std::list<CRSInfo> getCRSInfoList() const;
+
     // non-standard
     PROJ_DLL static AuthorityFactoryNNPtr
     create(const DatabaseContextNNPtr &context,
diff --git a/m4/pkg.m4 b/m4/pkg.m4
new file mode 100644
index 00000000..452488c8
--- /dev/null
+++ b/m4/pkg.m4
@@ -0,0 +1,157 @@
+# pkg.m4 - Macros to locate and utilize pkg-config.            -*- Autoconf -*-
+# 
+# Copyright © 2004 Scott James Remnant <scott@netsplit.com>.
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+#
+# As a special exception to the GNU General Public License, if you
+# distribute this file as part of a program that contains a
+# configuration script generated by Autoconf, you may include it under
+# the same distribution terms that you use for the rest of that program.
+
+# PKG_PROG_PKG_CONFIG([MIN-VERSION])
+# ----------------------------------
+AC_DEFUN([PKG_PROG_PKG_CONFIG],
+[m4_pattern_forbid([^_?PKG_[A-Z_]+$])
+m4_pattern_allow([^PKG_CONFIG(_PATH)?$])
+AC_ARG_VAR([PKG_CONFIG], [path to pkg-config utility])dnl
+if test "x$ac_cv_env_PKG_CONFIG_set" != "xset"; then
+	AC_PATH_TOOL([PKG_CONFIG], [pkg-config])
+fi
+if test -n "$PKG_CONFIG"; then
+	_pkg_min_version=m4_default([$1], [0.9.0])
+	AC_MSG_CHECKING([pkg-config is at least version $_pkg_min_version])
+	if $PKG_CONFIG --atleast-pkgconfig-version $_pkg_min_version; then
+		AC_MSG_RESULT([yes])
+	else
+		AC_MSG_RESULT([no])
+		PKG_CONFIG=""
+	fi
+		
+fi[]dnl
+])# PKG_PROG_PKG_CONFIG
+
+# PKG_CHECK_EXISTS(MODULES, [ACTION-IF-FOUND], [ACTION-IF-NOT-FOUND])
+#
+# Check to see whether a particular set of modules exists.  Similar
+# to PKG_CHECK_MODULES(), but does not set variables or print errors.
+#
+#
+# Similar to PKG_CHECK_MODULES, make sure that the first instance of
+# this or PKG_CHECK_MODULES is called, or make sure to call
+# PKG_CHECK_EXISTS manually
+# --------------------------------------------------------------
+AC_DEFUN([PKG_CHECK_EXISTS],
+[AC_REQUIRE([PKG_PROG_PKG_CONFIG])dnl
+if test -n "$PKG_CONFIG" && \
+    AC_RUN_LOG([$PKG_CONFIG --exists --print-errors "$1"]); then
+  m4_ifval([$2], [$2], [:])
+m4_ifvaln([$3], [else
+  $3])dnl
+fi])
+
+
+# _PKG_CONFIG([VARIABLE], [COMMAND], [MODULES])
+# ---------------------------------------------
+m4_define([_PKG_CONFIG],
+[if test -n "$PKG_CONFIG"; then
+    if test -n "$$1"; then
+        pkg_cv_[]$1="$$1"
+    else
+        PKG_CHECK_EXISTS([$3],
+                         [pkg_cv_[]$1=`$PKG_CONFIG --[]$2 "$3" 2>/dev/null`],
+			 [pkg_failed=yes])
+    fi
+else
+	pkg_failed=untried
+fi[]dnl
+])# _PKG_CONFIG
+
+# _PKG_SHORT_ERRORS_SUPPORTED
+# -----------------------------
+AC_DEFUN([_PKG_SHORT_ERRORS_SUPPORTED],
+[AC_REQUIRE([PKG_PROG_PKG_CONFIG])
+if $PKG_CONFIG --atleast-pkgconfig-version 0.20; then
+        _pkg_short_errors_supported=yes
+else
+        _pkg_short_errors_supported=no
+fi[]dnl
+])# _PKG_SHORT_ERRORS_SUPPORTED
+
+
+# PKG_CHECK_MODULES(VARIABLE-PREFIX, MODULES, [ACTION-IF-FOUND],
+# [ACTION-IF-NOT-FOUND])
+#
+#
+# Note that if there is a possibility the first call to
+# PKG_CHECK_MODULES might not happen, you should be sure to include an
+# explicit call to PKG_PROG_PKG_CONFIG in your configure.ac
+#
+#
+# --------------------------------------------------------------
+AC_DEFUN([PKG_CHECK_MODULES],
+[AC_REQUIRE([PKG_PROG_PKG_CONFIG])dnl
+AC_ARG_VAR([$1][_CFLAGS], [C compiler flags for $1, overriding pkg-config])dnl
+AC_ARG_VAR([$1][_LIBS], [linker flags for $1, overriding pkg-config])dnl
+
+pkg_failed=no
+AC_MSG_CHECKING([for $1])
+
+_PKG_CONFIG([$1][_CFLAGS], [cflags], [$2])
+_PKG_CONFIG([$1][_LIBS], [libs], [$2])
+
+m4_define([_PKG_TEXT], [Alternatively, you may set the environment variables $1[]_CFLAGS
+and $1[]_LIBS to avoid the need to call pkg-config.
+See the pkg-config man page for more details.])
+
+if test $pkg_failed = yes; then
+        _PKG_SHORT_ERRORS_SUPPORTED
+        if test $_pkg_short_errors_supported = yes; then
+	        $1[]_PKG_ERRORS=`$PKG_CONFIG --short-errors --errors-to-stdout --print-errors "$2"`
+        else 
+	        $1[]_PKG_ERRORS=`$PKG_CONFIG --errors-to-stdout --print-errors "$2"`
+        fi
+	# Put the nasty error message in config.log where it belongs
+	echo "$$1[]_PKG_ERRORS" >&AS_MESSAGE_LOG_FD
+
+	ifelse([$4], , [AC_MSG_ERROR(dnl
+[Package requirements ($2) were not met:
+
+$$1_PKG_ERRORS
+
+Consider adjusting the PKG_CONFIG_PATH environment variable if you
+installed software in a non-standard prefix.
+
+_PKG_TEXT
+])],
+		[AC_MSG_RESULT([no])
+                $4])
+elif test $pkg_failed = untried; then
+	ifelse([$4], , [AC_MSG_FAILURE(dnl
+[The pkg-config script could not be found or is too old.  Make sure it
+is in your PATH or set the PKG_CONFIG environment variable to the full
+path to pkg-config.
+
+_PKG_TEXT
+
+To get pkg-config, see <http://pkg-config.freedesktop.org/>.])],
+		[$4])
+else
+	$1[]_CFLAGS=$pkg_cv_[]$1[]_CFLAGS
+	$1[]_LIBS=$pkg_cv_[]$1[]_LIBS
+        AC_MSG_RESULT([yes])
+	ifelse([$3], , :, [$3])
+fi[]dnl
+])# PKG_CHECK_MODULES
diff --git a/man/man1/nad2nad.1 b/man/man1/nad2nad.1
deleted file mode 100644
index 96db4049..00000000
--- a/man/man1/nad2nad.1
+++ /dev/null
@@ -1,195 +0,0 @@
-.\" release 5
-.\"
-.ad b
-.hy 1
-.TH NAD2NAD 1 "2018/02/15 Rel. 5.0.0"
-.SH NAME
-nad2nad \- North American Datum conversion filter
-.SH SYNOPSIS
-.B nad2nad
-[
-.B \-eEfihortwW
-[
-.I args
-] ] [
-.B +args
-]
-file[s]
-.SH DESCRIPTION
-Program
-.B nad2nad
-is a filter to convert data between North
-America Datum 1927 (NAD27) and North American Datum 1983.
-.B nad2nad
-can optionally process both State Plane Coordinate System (SPCS) and
-Universal Transverse Mercator (UTM) grid data as well
-as geographic data for both input and output.  This can also be 
-accomplished with the \fIcs2cs\fR program.
-.PP
-The following control parameters can appear in any order:
-.TP
-.BI \-[i|o] " keyword[,keyword]"
-The
-.B \-i
-and
-.B \-o
-option expect keyword arguments which define various characteristics and
-processing modes of the
-respective input data.
-Usage allows multiple arguments to be included with a \- operator
-when separated by commas.
-Datum conversion requires the data to be in geographic coordinates, but
-.B nad2nad
-will allow conversion of data to and from SPCS or UTM grid systems.
-The following are keywords and arguments recognized by both the
-.B \-i
-and
-.B \-o
-that will apply to respective input and output conversion of user
-data to internal geographic coordinates:
-.RS 1in
-.TP
-.B 27|83
-datum of data
-.TP
-.BI utm= n
-.B UTM
-coordinates in meters for zone
-.I n
-.TP
-.BI spcs= n
-for data in
-.B SPCS
-coordinates, where
-.I n
-is state zone number.
-.TP
-.B feet
-data units are in U.S. Surveyor's feet.
-This is allowed only when the
-.I spcs
-option been previously used.
-Default coordinates are in meters.
-.TP
-.B bin
-for data in binary form.
-.TP
-.B rev
-data in lat-lon order rather than default lon-lat order.
-.TP
-.BI hp= ss
-use high precision conversion zone
-.I ss.
-Certain States have ancillary correction tables to further refine
-the basic conus table.
-.I Ss
-key and States are:
-.RS .5in
-FL Florida
-.br
-MD Maryland
-.br
-TN Tennessee
-.br
-WI Wisconsin
-.br
-WO Washington, Oregon and northern part of California.
-.RE
-.RE
-.TP
-.BI \-t "a"
-.I A
-specifies a character employed as the first character to denote
-a control line to be passed through without processing.
-This option applicable to ascii input only.
-(# is the default value).
-.TP
-.BI \-e " string"
-.I String
-is an arbitrary string to be output if an error is detected during
-data transformations.
-The default value is: *\et*.
-Note that if the
-.B "\-o bin"
-option is employed, an error is output as HUGE_VAL for both values.
-.TP
-.BI \-r " region"
-specifies which regional conversion table to employ which are identified
-by the following:
-.RS .5in
-conus \- conterminous 48 States
-.br
-alaska \- State of Alaska
-.br
-hawaii \- State of Hawaii
-.br
-prvi \- Puerto Rico and Virgin Islands
-.br
-stgeorge \- St. George Is, Alaska
-.br
-stpaul \- St. Paul Is, Alaska
-.br
-stlrnc \- St. Lawrence Is, Alaska
-.TP
-.B \-E
-Input coordinates are echoed to output before output values.
-.RE
-.TP
-.BI \-f " format"
-.I Format
-is a
-.I printf
-format string to control the form of the output values.
-For inverse projections, the output will be in degrees when this option
-is employed.
-If a format is specified for inverse projection the
-output data will be in decimal degrees.
-The default format is \(``%.2f\('' for forward projection and DMS
-for inverse.
-.TP
-.BI \-[w|W] n
-.I N
-is the number of significant fractional digits to employ for
-seconds output (when the option is not specified,
-.B \-w3
-is assumed).
-When
-.B \-W
-is employed the fields will be constant width and with leading zeroes.
-.PP
-One or more
-.I files
-(processed in left to right order)
-specify the source of data to be transformed.
-A \- will specify the location of processing standard input.
-If no files are specified, the input is assumed to be from
-.I stdin.
-For ASCII input data the two data values must be in the
-first two white space separated fields and
-when both input and output are ASCII all trailing portions
-of the input line are appended to the output line.
-.PP
-Input geographic data
-(longitude and latitude) must be in DMS format when neither
-.I utm
-nor
-.I spcs
-is specified, otherwise in meters or feet (\fIfeet\fR option used).
-Input data fields must be separated by white space and not have
-imbedded white space.
-.PP
-Output data will be in tab separated fields of DMS or grid
-coordinates in meters or feet.
-.PP
-Any data after the two input values are echoed after the two
-output data values.
-.SH SEE ALSO
-.B cs2cs(1U), proj(1U),
-.br
-.I "Cartographic Projection Procedures for the UNIX Environment\(emA User's Manual,"
-(Evenden, 1990, Open-file report 90\-284).
-.SH BUGS
-A list of known bugs can found at https://github.com/OSGeo/proj.4/issues
-where new bug reports can be submitted too.
-.SH HOME PAGE
-https://proj4.org/
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index de8ca3e5..122227bf 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -7,7 +7,6 @@ option(BUILD_CCT     "Build cct (coordinate conversion and transformation tool)"
 option(BUILD_CS2CS   "Build cs2cs (coordinate systems to coordinate systems translation tool)" ON)
 option(BUILD_GEOD    "Build geod (computation of geodesic lines)" ON)
 option(BUILD_GIE     "Build gie (geospatial integrity investigation environment - a PROJ.4 test tool)" ON)
-option(BUILD_NAD2BIN "Build nad2bin (format conversion tool)" ON)
 option(BUILD_PROJ    "Build proj (cartographic projection tool : latlong <-> projected coordinates)" ON)
 option(BUILD_PROJINFO "Build projinfo (SRS and coordinate operation metadata/query tool)" ON)
 
@@ -50,11 +49,6 @@ if(BUILD_GEOD)
   set(BIN_TARGETS ${BIN_TARGETS} geod)
 endif(BUILD_GEOD)
 
-if(BUILD_NAD2BIN)
-  include(bin_nad2bin.cmake)
-  set(BIN_TARGETS ${BIN_TARGETS} nad2bin)
-endif(BUILD_NAD2BIN)
-
 if(BUILD_PROJ)
   include(bin_proj.cmake)
   set(BIN_TARGETS ${BIN_TARGETS} binproj)
diff --git a/src/Makefile.am b/src/Makefile.am
index fabea4bf..698baf97 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,20 +1,20 @@
 AM_CFLAGS = @C_WFLAGS@
 
-bin_PROGRAMS =	proj nad2bin geod cs2cs gie cct projinfo
+bin_PROGRAMS =	proj geod cs2cs gie cct projinfo
 EXTRA_PROGRAMS = multistresstest test228
 
 TESTS = geodtest
 check_PROGRAMS = geodtest
 
 AM_CPPFLAGS =	-DPROJ_LIB=\"$(pkgdatadir)\" \
-		-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@ -I$(top_srcdir)/include @SQLITE3_FLAGS@
+		-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@ -I$(top_srcdir)/include @SQLITE3_CFLAGS@
 AM_CXXFLAGS =    @CXX_WFLAGS@ @FLTO_FLAG@ -DPROJ_COMPILATION
 
 include_HEADERS = proj.h proj_experimental.h proj_constants.h proj_api.h geodesic.h \
 	org_proj4_PJ.h proj_symbol_rename.h
 
 EXTRA_DIST = bin_cct.cmake bin_gie.cmake bin_cs2cs.cmake \
-	bin_geod.cmake bin_nad2bin.cmake bin_proj.cmake bin_projinfo.cmake \
+	bin_geod.cmake bin_proj.cmake bin_projinfo.cmake \
 	lib_proj.cmake CMakeLists.txt bin_geodtest.cmake tests/geodtest.cpp \
 	wkt1_grammar.y wkt2_grammar.y apps/emess.h
 
@@ -22,7 +22,6 @@ proj_SOURCES = apps/proj.cpp apps/emess.cpp
 projinfo_SOURCES = apps/projinfo.cpp
 cs2cs_SOURCES = apps/cs2cs.cpp apps/emess.cpp
 cct_SOURCES = apps/cct.cpp apps/proj_strtod.cpp apps/proj_strtod.h apps/optargpm.h
-nad2bin_SOURCES = apps/nad2bin.cpp
 geod_SOURCES = apps/geod.cpp apps/geod_set.cpp apps/geod_interface.cpp apps/geod_interface.h apps/emess.cpp
 
 gie_SOURCES = apps/gie.cpp apps/proj_strtod.cpp apps/proj_strtod.h apps/optargpm.h
@@ -35,7 +34,6 @@ cs2cs_LDADD = libproj.la
 geod_LDADD = libproj.la
 proj_LDADD = libproj.la
 projinfo_LDADD = libproj.la
-nad2bin_LDADD = libproj.la
 
 gie_LDADD = libproj.la
 multistresstest_LDADD = libproj.la @THREAD_LIB@
@@ -45,7 +43,7 @@ geodtest_LDADD = libproj.la
 lib_LTLIBRARIES = libproj.la
 
 libproj_la_LDFLAGS = -no-undefined -version-info 14:1:1
-libproj_la_LIBADD = @SQLITE3_LDFLAGS@
+libproj_la_LIBADD = @SQLITE3_LIBS@
 
 libproj_la_SOURCES = \
 	pj_list.h proj_internal.h proj_math.h \
@@ -124,7 +122,6 @@ libproj_la_SOURCES = \
 	projections/wag7.cpp \
 	projections/lcca.cpp \
 	projections/geos.cpp \
-	projections/etmerc.cpp \
 	projections/boggs.cpp \
 	projections/collg.cpp \
 	projections/comill.cpp \
@@ -186,15 +183,15 @@ libproj_la_SOURCES = \
 	transformations/molodensky.cpp \
 	transformations/vgridshift.cpp \
 	\
-	aasincos.cpp adjlon.cpp bch2bps.cpp bchgen.cpp \
-	biveval.cpp dmstor.cpp mk_cheby.cpp auth.cpp \
+	aasincos.cpp adjlon.cpp \
+	dmstor.cpp auth.cpp \
 	deriv.cpp ell_set.cpp ellps.cpp errno.cpp \
 	factors.cpp fwd.cpp init.cpp inv.cpp \
 	list.cpp malloc.cpp mlfn.cpp msfn.cpp proj_mdist.cpp \
 	open_lib.cpp param.cpp phi2.cpp pr_list.cpp \
 	qsfn.cpp strerrno.cpp \
 	tsfn.cpp units.cpp ctx.cpp log.cpp zpoly1.cpp rtodms.cpp \
-	vector1.cpp release.cpp gauss.cpp \
+	release.cpp gauss.cpp \
 	fileapi.cpp \
 	\
 	gc_reader.cpp gridcatalog.cpp \
diff --git a/src/apps/cct.cpp b/src/apps/cct.cpp
index 4deefba6..34bf0777 100644
--- a/src/apps/cct.cpp
+++ b/src/apps/cct.cpp
@@ -410,6 +410,8 @@ int main(int argc, char **argv) {
             );
     }
 
+    proj_destroy(P);
+
     if (stdout != fout)
         fclose (fout);
     free (o);
diff --git a/src/apps/gie.cpp b/src/apps/gie.cpp
index e912a076..c3622d52 100644
--- a/src/apps/gie.cpp
+++ b/src/apps/gie.cpp
@@ -1058,7 +1058,7 @@ static const struct errno_vs_err_const lookup[] = {
     {"pjd_err_lat_0_or_alpha_eq_90"     ,  -33},
     {"pjd_err_ellipsoid_use_required"   ,  -34},
     {"pjd_err_invalid_utm_zone"         ,  -35},
-    {"pjd_err_tcheby_val_out_of_range"  ,  -36},
+    {""                                 ,  -36}, /* no longer used */
     {"pjd_err_failed_to_find_proj"      ,  -37},
     {"pjd_err_failed_to_load_grid"      ,  -38},
     {"pjd_err_invalid_m_or_n"           ,  -39},
@@ -1081,6 +1081,7 @@ static const struct errno_vs_err_const lookup[] = {
     {"pjd_err_ellipsoidal_unsupported"  ,  -56},
     {"pjd_err_too_many_inits"           ,  -57},
     {"pjd_err_invalid_arg"              ,  -58},
+    {"pjd_err_inconsistent_unit"        ,  -59},
     {"pjd_err_dont_skip"                ,  5555},
     {"pjd_err_unknown"                  ,  9999},
     {"pjd_err_enomem"                   ,  ENOMEM},
@@ -1140,7 +1141,8 @@ static int errno_from_err_const (const char *err_const) {
 
     /* First try to find a match excluding the PJD_ERR_ prefix */
     for (i = 0;  i < n;  i++) {
-        if (0==strncmp (lookup[i].the_err_const + 8, err_const, len))
+        if (strlen(lookup[i].the_err_const) > 8 &&
+            0==strncmp (lookup[i].the_err_const + 8, err_const, len))
             return lookup[i].the_errno;
     }
 
diff --git a/src/apps/nad2bin.cpp b/src/apps/nad2bin.cpp
deleted file mode 100644
index a684b087..00000000
--- a/src/apps/nad2bin.cpp
+++ /dev/null
@@ -1,382 +0,0 @@
-/* Convert bivariate ASCII NAD27 to NAD83 tables to NTv2 binary structure */
-#include <stdio.h>
-#include <stdlib.h>
-
-#define PJ_LIB__
-#include "proj_internal.h"
-#include "proj_internal.h"
-#define U_SEC_TO_RAD 4.848136811095359935899141023e-12
-
-/************************************************************************/
-/*                             swap_words()                             */
-/*                                                                      */
-/*      Convert the byte order of the given word(s) in place.           */
-/************************************************************************/
-
-static const int  byte_order_test = 1;
-#define IS_LSB	(((const unsigned char *) (&byte_order_test))[0] == 1)
-
-static void swap_words( void *data_in, int word_size, int word_count )
-
-{
-    int	word;
-    unsigned char *data = (unsigned char *) data_in;
-
-    for( word = 0; word < word_count; word++ )
-    {
-        int	i;
-        
-        for( i = 0; i < word_size/2; i++ )
-        {
-            unsigned char	t;
-            
-            t = data[i];
-            data[i] = data[word_size-i-1];
-            data[word_size-i-1] = t;
-        }
-        
-        data += word_size;
-    }
-}
-
-/************************************************************************/
-/*                               Usage()                                */
-/************************************************************************/
-
-static void Usage()
-{
-    fprintf(stderr,
-            "usage: nad2bin [-f ctable/ctable2/ntv2] binary_output < ascii_source\n" );
-    exit(1);
-}
-
-/************************************************************************/
-/*                                main()                                */
-/************************************************************************/
-int main(int argc, char **argv) {
-    struct CTABLE ct;
-    FLP *p, t;
-    size_t tsize;
-    int i, j, ichk;
-    long lam, laml, phi, phil;
-    FILE *fp;
-
-    const char *output_file = nullptr;
-
-    const char *format   = "ctable2";
-    const char *GS_TYPE  = "SECONDS";
-    const char *VERSION  = "";
-    const char *SYSTEM_F = "NAD27";
-    const char *SYSTEM_T = "NAD83";
-    const char *SUB_NAME = "";
-    const char *CREATED  = "";
-    const char *UPDATED  = "";
-
-/* ==================================================================== */
-/*      Process arguments.                                              */
-/* ==================================================================== */
-    for( i = 1; i < argc; i++ )
-    {
-        if( i < argc-1 && strcmp(argv[i],"-f") == 0 ) 
-        {
-            format = argv[++i];
-        }
-        else if( output_file == nullptr )
-        {
-            output_file = argv[i];
-        }
-        else
-            Usage();
-    }
-
-    if( output_file == nullptr )
-        Usage();
-
-    fprintf( stdout, "Output Binary File Format: %s\n", format );
-
-/* ==================================================================== */
-/*      Read the ASCII Table                                            */
-/* ==================================================================== */
-
-    memset(ct.id,0,MAX_TAB_ID);
-    if ( nullptr == fgets(ct.id, MAX_TAB_ID, stdin) ) {
-        perror("fgets");
-        exit(1);
-    }
-    /* cppcheck-suppress invalidscanf */
-    if ( EOF == scanf("%d %d %*d %lf %lf %lf %lf", &ct.lim.lam, &ct.lim.phi,
-          &ct.ll.lam, &ct.del.lam, &ct.ll.phi, &ct.del.phi) ) {
-        perror("scanf");
-        exit(1);
-    }
-    if (!(ct.cvs = (FLP *)malloc(tsize = ct.lim.lam * ct.lim.phi *
-                                 sizeof(FLP)))) {
-        perror("mem. alloc");
-        exit(1);
-    }
-    ct.ll.lam *= DEG_TO_RAD;
-    ct.ll.phi *= DEG_TO_RAD;
-    ct.del.lam *= DEG_TO_RAD;
-    ct.del.phi *= DEG_TO_RAD;
-    /* load table */
-    p = ct.cvs;
-    for (i = 0; i < ct.lim.phi; ++i) {
-        /* cppcheck-suppress invalidscanf */
-        if ( EOF == scanf("%d:%ld %ld", &ichk, &laml, &phil) ) {
-            perror("scanf on row");
-            exit(1);
-        }
-        if (ichk != i) {
-            fprintf(stderr,"format check on row\n");
-            exit(1);
-        }
-        t.lam = (float) (laml * U_SEC_TO_RAD);
-        t.phi = (float) (phil * U_SEC_TO_RAD);
-        *p++ = t;
-        for (j = 1; j < ct.lim.lam; ++j) {
-            /* cppcheck-suppress invalidscanf */
-            if ( EOF == scanf("%ld %ld", &lam, &phi) ) {
-                perror("scanf on column");
-                exit(1);
-            }
-            t.lam = (float) ((laml += lam) * U_SEC_TO_RAD);
-            t.phi = (float) ((phil += phi) * U_SEC_TO_RAD);
-            *p++ = t;
-        }
-    }
-    if (feof(stdin)) {
-        fprintf(stderr, "premature EOF\n");
-        exit(1);
-    }
-
-/* ==================================================================== */
-/*      Write out the old ctable format - this is machine and byte      */
-/*      order specific.                                                 */
-/* ==================================================================== */
-    if( strcmp(format,"ctable") == 0 ) 
-    {
-	if (!(fp = fopen(output_file, "wb"))) {
-            perror(output_file);
-            exit(2);
-	}
-	if (fwrite(&ct, sizeof(ct), 1, fp) != 1 ||
-            fwrite(ct.cvs, tsize, 1, fp) != 1) {
-            fprintf(stderr, "output failure\n");
-            exit(2);
-	}
-        fclose( fp );
-	exit(0); /* normal completion */
-    }
-
-/* ==================================================================== */
-/*      Write out the old ctable format - this is machine and byte      */
-/*      order specific.                                                 */
-/* ==================================================================== */
-    if( strcmp(format,"ctable2") == 0 ) 
-    {
-        char header[160];
-
-	if (!(fp = fopen(output_file, "wb"))) {
-            perror(output_file);
-            exit(2);
-	}
-
-	/* cppcheck-suppress sizeofCalculation */
-        STATIC_ASSERT( MAX_TAB_ID == 80 );
-        /* cppcheck-suppress sizeofCalculation */
-        STATIC_ASSERT( sizeof(pj_int32) == 4 ); /* for ct.lim.lam/phi */
-
-        memset( header, 0, sizeof(header) );
-
-        memcpy( header +   0, "CTABLE V2.0     ", 16 );
-        memcpy( header +  16, ct.id, 80 );
-        memcpy( header +  96, &ct.ll.lam, 8 );
-        memcpy( header + 104, &ct.ll.phi, 8 );
-        memcpy( header + 112, &ct.del.lam, 8 );
-        memcpy( header + 120, &ct.del.phi, 8 );
-        memcpy( header + 128, &ct.lim.lam, 4 );
-        memcpy( header + 132, &ct.lim.phi, 4 );
-
-        /* force into LSB format */
-        if( !IS_LSB ) 
-        {
-            swap_words( header +  96, 8, 4 );
-            swap_words( header + 128, 4, 2 );
-            swap_words( ct.cvs, 4, ct.lim.lam * 2 * ct.lim.phi );
-        }
-
-        if( fwrite( header, sizeof(header), 1, fp ) != 1 ) {
-            perror( "fwrite" );
-            exit( 2 );
-        }
-
-	if (fwrite(ct.cvs, tsize, 1, fp) != 1) {
-            perror( "fwrite" );
-            exit(2);
-	}
-
-        fclose( fp );
-	exit(0); /* normal completion */
-    }
-
-/* ==================================================================== */
-/*      Write out the NTv2 format grid shift file.                      */
-/* ==================================================================== */
-    if( strcmp(format,"ntv2") == 0 ) 
-    {
-        if (!(fp = fopen(output_file, "wb"))) 
-        {
-            perror(output_file);
-            exit(2);
-        }
-        
-/* -------------------------------------------------------------------- */
-/*      Write the file header.                                          */
-/* -------------------------------------------------------------------- */
-        {    
-            char achHeader[11*16];
-
-            memset( achHeader, 0, sizeof(achHeader) );
-        
-            memcpy( achHeader +  0*16, "NUM_OREC", 8 );
-            achHeader[ 0*16 + 8] = 0xb;
-
-            memcpy( achHeader +  1*16, "NUM_SREC", 8 );
-            achHeader[ 1*16 + 8] = 0xb;
-
-            memcpy( achHeader +  2*16, "NUM_FILE", 8 );
-            achHeader[ 2*16 + 8] = 0x1;
-
-            memcpy( achHeader +  3*16, "GS_TYPE         ", 16 );
-            memcpy( achHeader +  3*16+8, GS_TYPE, MIN(16,strlen(GS_TYPE)) );
-
-            memcpy( achHeader +  4*16, "VERSION         ", 16 );
-            memcpy( achHeader +  4*16+8, VERSION, MIN(16,strlen(VERSION)) );
-
-            memcpy( achHeader +  5*16, "SYSTEM_F        ", 16 );
-            memcpy( achHeader +  5*16+8, SYSTEM_F, MIN(16,strlen(SYSTEM_F)) );
-
-            memcpy( achHeader +  6*16, "SYSTEM_T        ", 16 );
-            memcpy( achHeader +  6*16+8, SYSTEM_T, MIN(16,strlen(SYSTEM_T)) );
-
-            memcpy( achHeader +  7*16, "MAJOR_F ", 8);
-            memcpy( achHeader +  8*16, "MINOR_F ", 8 );
-            memcpy( achHeader +  9*16, "MAJOR_T ", 8 );
-            memcpy( achHeader + 10*16, "MINOR_T ", 8 );
-
-            fwrite( achHeader, 1, sizeof(achHeader), fp );
-        }
-        
-/* -------------------------------------------------------------------- */
-/*      Write the grid header.                                          */
-/* -------------------------------------------------------------------- */
-        {
-            unsigned char achHeader[11*16];
-            double dfValue;
-            pj_int32 nGSCount = ct.lim.lam * ct.lim.phi;
-            PJ_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;
-
-            /* cppcheck-suppress sizeofCalculation */
-            STATIC_ASSERT( sizeof(nGSCount) == 4 );
-
-            memset( achHeader, 0, sizeof(achHeader) );
-
-            memcpy( achHeader +  0*16, "SUB_NAME        ", 16 );
-            memcpy( achHeader +  0*16+8, SUB_NAME, MIN(16,strlen(SUB_NAME)) );
-    
-            memcpy( achHeader +  1*16, "PARENT          ", 16 );
-            memcpy( achHeader +  1*16+8, "NONE", MIN(16,strlen("NONE")) );
-    
-            memcpy( achHeader +  2*16, "CREATED         ", 16 );
-            memcpy( achHeader +  2*16+8, CREATED, MIN(16,strlen(CREATED)) );
-    
-            memcpy( achHeader +  3*16, "UPDATED         ", 16 );
-            memcpy( achHeader +  3*16+8, UPDATED, MIN(16,strlen(UPDATED)) );
-
-            memcpy( achHeader +  4*16, "S_LAT   ", 8 );
-            dfValue = ct.ll.phi * 3600.0 / DEG_TO_RAD;
-            memcpy( achHeader +  4*16 + 8, &dfValue, 8 );
-
-            memcpy( achHeader +  5*16, "N_LAT   ", 8 );
-            dfValue = ur.phi * 3600.0 / DEG_TO_RAD;
-            memcpy( achHeader +  5*16 + 8, &dfValue, 8 );
-
-            memcpy( achHeader +  6*16, "E_LONG  ", 8 );
-            dfValue = -1 * ur.lam * 3600.0 / DEG_TO_RAD;
-            memcpy( achHeader +  6*16 + 8, &dfValue, 8 );
-
-            memcpy( achHeader +  7*16, "W_LONG  ", 8 );
-            dfValue = -1 * ct.ll.lam * 3600.0 / DEG_TO_RAD;
-            memcpy( achHeader +  7*16 + 8, &dfValue, 8 );
-
-            memcpy( achHeader +  8*16, "LAT_INC ", 8 );
-            dfValue = ct.del.phi * 3600.0 / DEG_TO_RAD;
-            memcpy( achHeader +  8*16 + 8, &dfValue, 8 );
-    
-            memcpy( achHeader +  9*16, "LONG_INC", 8 );
-            dfValue = ct.del.lam * 3600.0 / DEG_TO_RAD;
-            memcpy( achHeader +  9*16 + 8, &dfValue, 8 );
-    
-            memcpy( achHeader + 10*16, "GS_COUNT", 8 );
-            memcpy( achHeader + 10*16+8, &nGSCount, 4 );
-    
-            if( !IS_LSB ) 
-            {
-                swap_words( achHeader +  4*16 + 8, 8, 1 );
-                swap_words( achHeader +  5*16 + 8, 8, 1 );
-                swap_words( achHeader +  6*16 + 8, 8, 1 );
-                swap_words( achHeader +  7*16 + 8, 8, 1 );
-                swap_words( achHeader +  8*16 + 8, 8, 1 );
-                swap_words( achHeader +  9*16 + 8, 8, 1 );
-                swap_words( achHeader + 10*16 + 8, 4, 1 );
-            }
-
-            fwrite( achHeader, 1, sizeof(achHeader), fp );
-        }
-
-/* -------------------------------------------------------------------- */
-/*      Write the actual grid cells.                                    */
-/* -------------------------------------------------------------------- */
-        {
-            float *row_buf;
-            int row;
-
-            row_buf = (float *) pj_malloc(ct.lim.lam * sizeof(float) * 4);
-            memset( row_buf, 0, sizeof(float)*4 );
-
-            for( row = 0; row < ct.lim.phi; row++ )
-            {
-                for( i = 0; i < ct.lim.lam; i++ )
-                {
-                    FLP *cvs = ct.cvs + (row) * ct.lim.lam
-                        + (ct.lim.lam - i - 1);
-
-                    /* convert radians to seconds */
-                    row_buf[i*4+0] = (float) (cvs->phi * (3600.0 / (M_PI/180.0)));
-                    row_buf[i*4+1] = (float) (cvs->lam * (3600.0 / (M_PI/180.0)));
-
-                    /* We leave the accuracy values as zero */
-                }
-
-                if( !IS_LSB )
-                    swap_words( row_buf, 4, ct.lim.lam * 4 );
-
-                if( fwrite( row_buf, sizeof(float), ct.lim.lam*4, fp ) 
-                    != (size_t)( 4 * ct.lim.lam ) )
-                {
-                    perror( "write()" );
-                    exit( 2 );
-                }
-            }
-        }
-
-        fclose( fp );
-        exit(0); /* normal completion */
-    }
-
-    fprintf( stderr, "Unsupported format, nothing written.\n" );
-    exit( 3 );
-}
diff --git a/src/bch2bps.cpp b/src/bch2bps.cpp
deleted file mode 100644
index 9346457c..00000000
--- a/src/bch2bps.cpp
+++ /dev/null
@@ -1,154 +0,0 @@
-/* convert bivariate w Chebyshev series to w Power series */
-
-#include "proj.h"
-#include "proj_internal.h"
-/* basic support procedures */
-	static void /* clear vector to zero */
-clear(PJ_UV *p, int n) { static const PJ_UV c = {0., 0.}; while (n--) *p++ = c; }
-	static void /* clear matrix rows to zero */
-bclear(PJ_UV **p, int n, int m) { while (n--) clear(*p++, m); }
-	static void /* move vector */
-bmove(PJ_UV *a, PJ_UV *b, int n) { while (n--) *a++ = *b++; }
-	static void /* a <- m * b - c */
-submop(PJ_UV *a, double m, PJ_UV *b, PJ_UV *c, int n) {
-	while (n--) {
-		a->u = m * b->u - c->u;
-		a++->v = m * b++->v - c++->v;
-	}
-}
-	static void /* a <- b - c */
-subop(PJ_UV *a, PJ_UV *b, PJ_UV *c, int n) {
-	while (n--) {
-		a->u = b->u - c->u;
-		a++->v = b++->v - c++->v;
-	}
-}
-	static void /* multiply vector a by scalar m */
-dmult(PJ_UV *a, double m, int n) { while(n--) { a->u *= m; a->v *= m; ++a; } }
-	static void /* row adjust a[] <- a[] - m * b[] */
-dadd(PJ_UV *a, PJ_UV *b, double m, int n) {
-	while(n--) {
-		a->u -= m * b->u;
-		a++->v -= m * b++->v;
-	}
-}
-	static int /* convert row to power series */
-rows(PJ_UV *c, PJ_UV *d, int n) {
-	PJ_UV sv, *dd;
-	int j, k;
-
-	dd = (PJ_UV *)vector1(n-1, sizeof(PJ_UV));
-	if (!dd)
-		return 0;
-	sv.u = sv.v = 0.;
-	for (j = 0; j < n; ++j) d[j] = dd[j] = sv;
-	d[0] = c[n-1];
-	for (j = n-2; j >= 1; --j) {
-		for (k = n-j; k >= 1; --k) {
-			sv = d[k];
-			d[k].u = 2. * d[k-1].u - dd[k].u;
-			d[k].v = 2. * d[k-1].v - dd[k].v;
-			dd[k] = sv;
-		}
-		sv = d[0];
-		d[0].u = -dd[0].u + c[j].u;
-		d[0].v = -dd[0].v + c[j].v;
-		dd[0] = sv;
-	}
-	for (j = n-1; j >= 1; --j) {
-		d[j].u = d[j-1].u - dd[j].u;
-		d[j].v = d[j-1].v - dd[j].v;
-	}
-	d[0].u = -dd[0].u + .5 * c[0].u;
-	d[0].v = -dd[0].v + .5 * c[0].v;
-	pj_dalloc(dd);
-	return 1;
-}
-	static int /* convert columns to power series */
-cols(PJ_UV **c, PJ_UV **d, int nu, int nv) {
-	PJ_UV *sv, **dd;
-	int j, k;
-
-	dd = (PJ_UV **)vector2(nu, nv, sizeof(PJ_UV));
-	if (!dd)
-		return 0;
-	sv = (PJ_UV *)vector1(nv, sizeof(PJ_UV));
-	if (!sv) {
-		freev2((void **)dd, nu);
-		return 0;
-	}
-	bclear(d, nu, nv);
-	bclear(dd, nu, nv);
-	bmove(d[0], c[nu-1], nv);
-	for (j = nu-2; j >= 1; --j) {
-		for (k = nu-j; k >= 1; --k) {
-			bmove(sv, d[k], nv);
-			submop(d[k], 2., d[k-1], dd[k], nv);
-			bmove(dd[k], sv, nv);
-		}
-		bmove(sv, d[0], nv);
-		subop(d[0], c[j], dd[0], nv);
-		bmove(dd[0], sv, nv);
-	}
-	for (j = nu-1; j >= 1; --j)
-		subop(d[j], d[j-1], dd[j], nv);
-	submop(d[0], .5, c[0], dd[0], nv);
-	freev2((void **) dd, nu);
-	pj_dalloc(sv);
-	return 1;
-}
-	static void /* row adjust for range -1 to 1 to a to b */
-rowshft(double a, double b, PJ_UV *d, int n) {
-	int k, j;
-	double fac, cnst;
-
-	cnst = 2. / (b - a);
-	fac = cnst;
-	for (j = 1; j < n; ++j) {
-		d[j].u *= fac;
-		d[j].v *= fac;
-		fac *= cnst;
-	}
-	cnst = .5 * (a + b);
-	for (j = 0; j <= n-2; ++j)
-		for (k = n - 2; k >= j; --k) {
-			d[k].u -= cnst * d[k+1].u;
-			d[k].v -= cnst * d[k+1].v;
-		}
-}
-	static void /* column adjust for range -1 to 1 to a to b */
-colshft(double a, double b, PJ_UV **d, int n, int m) {
-	int k, j;
-	double fac, cnst;
-
-	cnst = 2. / (b - a);
-	fac = cnst;
-	for (j = 1; j < n; ++j) {
-		dmult(d[j], fac, m);
-		fac *= cnst;
-	}
-	cnst = .5 * (a + b);
-	for (j = 0; j <= n-2; ++j)
-		for (k = n - 2; k >= j; --k)
-			dadd(d[k], d[k+1], cnst, m);
-}
-	int /* entry point */
-bch2bps(PJ_UV a, PJ_UV b, PJ_UV **c, int nu, int nv) {
-	PJ_UV **d;
-	int i;
-
-	if (nu < 1 || nv < 1 || !(d = (PJ_UV **)vector2(nu, nv, sizeof(PJ_UV))))
-		return 0;
-	/* do rows to power series */
-	for (i = 0; i < nu; ++i) {
-		if (!rows(c[i], d[i], nv))
-			return 0;
-		rowshft(a.v, b.v, d[i], nv);
-	}
-	/* do columns to power series */
-	if (!cols(d, c, nu, nv))
-		return 0;
-	colshft(a.u, b.u, c, nu, nv);
-	freev2((void **) d, nu);
-	return 1;
-}
diff --git a/src/bchgen.cpp b/src/bchgen.cpp
deleted file mode 100644
index 9677b6f2..00000000
--- a/src/bchgen.cpp
+++ /dev/null
@@ -1,59 +0,0 @@
-/* generate double bivariate Chebychev polynomial */
-#include "proj.h"
-#include "proj_internal.h"
-	int
-bchgen(PJ_UV a, PJ_UV b, int nu, int nv, PJ_UV **f, PJ_UV(*func)(PJ_UV)) {
-	int i, j, k;
-	PJ_UV arg, *t, bma, bpa, *c;
-	double d, fac;
-
-	bma.u = 0.5 * (b.u - a.u); bma.v = 0.5 * (b.v - a.v);
-	bpa.u = 0.5 * (b.u + a.u); bpa.v = 0.5 * (b.v + a.v);
-	for ( i = 0; i < nu; ++i) {
-		arg.u = cos(M_PI * (i + 0.5) / nu) * bma.u + bpa.u;
-		for ( j = 0; j < nv; ++j) {
-			arg.v = cos(M_PI * (j + 0.5) / nv) * bma.v + bpa.v;
-			f[i][j] = (*func)(arg);
-			if ((f[i][j]).u == HUGE_VAL)
-				return(1);
-		}
-	}
-	if (!(c = (PJ_UV *) vector1(nu, sizeof(PJ_UV)))) return 1;
-	fac = 2. / nu;
-	for ( j = 0; j < nv ; ++j) {
-		for ( i = 0; i < nu; ++i) {
-			arg.u = arg.v = 0.;
-			for (k = 0; k < nu; ++k) {
-				d = cos(M_PI * i * (k + .5) / nu);
-				arg.u += f[k][j].u * d;
-				arg.v += f[k][j].v * d;
-			}
-			arg.u *= fac;
-			arg.v *= fac;
-			c[i] = arg;
-		}
-		for (i = 0; i < nu; ++i)
-			f[i][j] = c[i];
-	}
-	pj_dalloc(c);
-	if (!(c = (PJ_UV*) vector1(nv, sizeof(PJ_UV)))) return 1;
-	fac = 2. / nv;
-	for ( i = 0; i < nu; ++i) {
-		t = f[i];
-		for (j = 0; j < nv; ++j) {
-			arg.u = arg.v = 0.;
-			for (k = 0; k < nv; ++k) {
-				d = cos(M_PI * j * (k + .5) / nv);
-				arg.u += t[k].u * d;
-				arg.v += t[k].v * d;
-			}
-			arg.u *= fac;
-			arg.v *= fac;
-			c[j] = arg;
-		}
-		f[i] = c;
-		c = t;
-	}
-	pj_dalloc(c);
-	return(0);
-}
diff --git a/src/bin_nad2bin.cmake b/src/bin_nad2bin.cmake
deleted file mode 100644
index 271aac93..00000000
--- a/src/bin_nad2bin.cmake
+++ /dev/null
@@ -1,17 +0,0 @@
-if(WIN32 AND BUILD_LIBPROJ_SHARED)
-    message(warning " nad2nad can't be build with a DLL proj library you need a static lib")
-endif(WIN32 AND BUILD_LIBPROJ_SHARED)
-
-
-set(NAD2BIN_SRC apps/nad2bin.cpp)
-source_group("Source Files\\Bin" FILES ${NAD2BIN_SRC})
-
-#Executable
-add_executable(nad2bin ${NAD2BIN_SRC})
-target_link_libraries(nad2bin ${PROJ_LIBRARIES})
-install(TARGETS nad2bin 
-        RUNTIME DESTINATION ${BINDIR})
-
-if(MSVC AND BUILD_LIBPROJ_SHARED)
-    target_compile_definitions(nad2bin PRIVATE PROJ_MSVC_DLL_IMPORT=1)
-endif()
diff --git a/src/biveval.cpp b/src/biveval.cpp
deleted file mode 100644
index 9ead2fb7..00000000
--- a/src/biveval.cpp
+++ /dev/null
@@ -1,88 +0,0 @@
-/* procedures for evaluating Tseries */
-#include "proj.h"
-#include "proj_internal.h"
-# define NEAR_ONE	1.00001
-static double ceval(struct PW_COEF *C, int n, PJ_UV w, PJ_UV w2) {
-    double d=0, dd=0, vd, vdd, tmp, *c;
-    int j;
-
-    for (C += n ; n-- ; --C ) {
-        if ( (j = C->m) != 0) {
-            vd = vdd = 0.;
-            for (c = C->c + --j; j ; --j ) {
-                vd = w2.v * (tmp = vd) - vdd + *c--;
-                vdd = tmp;
-            }
-            d = w2.u * (tmp = d) - dd + w.v * vd - vdd + 0.5 * *c;
-        } else
-            d = w2.u * (tmp = d) - dd;
-        dd = tmp;
-    }
-    if ( (j = C->m) != 0 ) {
-        vd = vdd = 0.;
-        for (c = C->c + --j; j ; --j ) {
-            vd = w2.v * (tmp = vd) - vdd + *c--;
-            vdd = tmp;
-        }
-        return (w.u * d - dd + 0.5 * ( w.v * vd - vdd + 0.5 * *c ));
-    } else
-        return (w.u * d - dd);
-}
-
-PJ_UV /* bivariate Chebyshev polynomial entry point */
-bcheval(PJ_UV in, Tseries *T) {
-    PJ_UV w2, w;
-    PJ_UV out;
-    /* scale to +-1 */
-    w.u = ( in.u + in.u - T->a.u ) * T->b.u;
-    w.v = ( in.v + in.v - T->a.v ) * T->b.v;
-    if (fabs(w.u) > NEAR_ONE || fabs(w.v) > NEAR_ONE) {
-        out.u = out.v = HUGE_VAL;
-        pj_errno = PJD_ERR_TCHEBY_VAL_OUT_OF_RANGE;
-    } else { /* double evaluation */
-        w2.u = w.u + w.u;
-        w2.v = w.v + w.v;
-        out.u = ceval(T->cu, T->mu, w, w2);
-        out.v = ceval(T->cv, T->mv, w, w2);
-    }
-    return out;
-}
-
-PJ_UV /* bivariate power polynomial entry point */
-bpseval(PJ_UV in, Tseries *T) {
-    PJ_UV out;
-    double *c, row;
-    int i, m;
-
-    out.u = out.v = 0.;
-    for (i = T->mu; i >= 0; --i) {
-        row = 0.;
-        if ((m = T->cu[i].m) != 0) {
-            c = T->cu[i].c + m;
-            while (m--)
-                row = *--c + in.v * row;
-        }
-        out.u = row + in.u * out.u;
-    }
-    for (i = T->mv; i >= 0; --i) {
-        row = 0.;
-        if ((m = T->cv[i].m) != 0) {
-            c = T->cv[i].c + m;
-            while (m--)
-                row = *--c + in.v * row;
-        }
-        out.v = row + in.u * out.v;
-    }
-    return out;
-}
-
-PJ_UV /* general entry point selecting evaluation mode */
-biveval(PJ_UV in, Tseries *T) {
-
-    if (T->power) {
-        return bpseval(in, T);
-    } else {
-        return bcheval(in, T);
-    }
-}
-
diff --git a/src/ell_set.cpp b/src/ell_set.cpp
index f4228be8..4c9fc892 100644
--- a/src/ell_set.cpp
+++ b/src/ell_set.cpp
@@ -88,7 +88,8 @@ int pj_ellipsoid (PJ *P) {
 
     /* Specifying R overrules everything */
     if (pj_get_param (P->params, "R")) {
-        ellps_size (P);
+        if (0 != ellps_size (P))
+            return 1;
         pj_calc_ellipsoid_params (P, P->a, 0);
         if (proj_errno (P))
             return 1;
diff --git a/src/init.cpp b/src/init.cpp
index 482e398d..2961bcca 100644
--- a/src/init.cpp
+++ b/src/init.cpp
@@ -764,20 +764,18 @@ pj_init_ctx_with_allow_init_epsg(projCtx ctx, int argc, char **argv, int allow_i
         s = units[i].to_meter;
     }
     if (s || (s = pj_param(ctx, start, "sto_meter").s)) {
-        double factor;
-        int ratio = 0;
-
-        /* ratio number? */
-        if (strlen (s) > 1 && s[0] == '1' && s[1]=='/') {
-            ratio = 1;
-            s += 2;
+        char* end_ptr = const_cast<char*>(s);
+        PIN->to_meter = pj_strtod(s, &end_ptr);
+        s = end_ptr;
+        if (*s == '/') { /* ratio number */
+            ++s;
+            double denom = pj_strtod(s, nullptr);
+            if (denom == 0.0)
+                return pj_default_destructor (PIN, PJD_ERR_UNIT_FACTOR_LESS_THAN_0);
+            PIN->to_meter /= denom;
         }
-
-        factor = pj_strtod(s, nullptr);
-        if ((factor <= 0.0) || (1/factor==0))
+        if (PIN->to_meter <= 0.0)
             return pj_default_destructor (PIN, PJD_ERR_UNIT_FACTOR_LESS_THAN_0);
-
-        PIN->to_meter = ratio?  1 / factor: factor;
         PIN->fr_meter = 1 / PIN->to_meter;
 
     } else
@@ -792,9 +790,16 @@ pj_init_ctx_with_allow_init_epsg(projCtx ctx, int argc, char **argv, int allow_i
         s = units[i].to_meter;
     }
     if (s || (s = pj_param(ctx, start, "svto_meter").s)) {
-        PIN->vto_meter = pj_strtod(s, nullptr);
-        if (*s == '/') /* ratio number */
-            PIN->vto_meter /= pj_strtod(++s, nullptr);
+        char* end_ptr = const_cast<char*>(s);
+        PIN->vto_meter = pj_strtod(s, &end_ptr);
+        s = end_ptr;
+        if (*s == '/') { /* ratio number */
+            ++s;
+            double denom = pj_strtod(s, nullptr);
+            if (denom == 0.0)
+                return pj_default_destructor (PIN, PJD_ERR_UNIT_FACTOR_LESS_THAN_0);
+            PIN->vto_meter /= denom;
+        }
         if (PIN->vto_meter <= 0.0)
             return pj_default_destructor (PIN, PJD_ERR_UNIT_FACTOR_LESS_THAN_0);
         PIN->vfr_meter = 1. / PIN->vto_meter;
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index de11f181..240d11b3 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -34,6 +34,7 @@
 #include <cstdarg>
 #include <cstring>
 #include <map>
+#include <new>
 #include <utility>
 #include <vector>
 
@@ -144,9 +145,11 @@ static PJ *pj_obj_create(PJ_CONTEXT *ctx, const IdentifiedObjectNNPtr &objIn) {
             // PROJ string.
         }
     }
-    auto pj = new PJ();
-    pj->descr = "ISO-19111 object";
-    pj->iso_obj = objIn;
+    auto pj = pj_new();
+    if (pj) {
+        pj->descr = "ISO-19111 object";
+        pj->iso_obj = objIn;
+    }
     return pj;
 }
 //! @endcond
@@ -1204,9 +1207,8 @@ const char *proj_as_wkt(PJ_CONTEXT *ctx, const PJ *obj, PJ_WKT_TYPE type,
  * @param type PROJ String version.
  * @param options NULL-terminated list of strings with "KEY=VALUE" format. or
  * NULL.
- * The currently recognized option is USE_ETMERC=YES to use
- * +proj=etmerc instead of +proj=tmerc (or USE_ETMERC=NO to disable implicit
- * use of etmerc by utm conversions)
+ * The currently recognized option is USE_APPROX_TMERC=YES to add the +approx
+ * flag to +proj=tmerc or +proj=utm
  * @return a string, or NULL in case of error.
  */
 const char *proj_as_proj_string(PJ_CONTEXT *ctx, const PJ *obj,
@@ -1240,10 +1242,8 @@ const char *proj_as_proj_string(PJ_CONTEXT *ctx, const PJ *obj,
     try {
         auto formatter = PROJStringFormatter::create(convention, dbContext);
         if (options != nullptr && options[0] != nullptr) {
-            if (ci_equal(options[0], "USE_ETMERC=YES")) {
-                formatter->setUseETMercForTMerc(true);
-            } else if (ci_equal(options[0], "USE_ETMERC=NO")) {
-                formatter->setUseETMercForTMerc(false);
+            if (ci_equal(options[0], "USE_APPROX_TMERC=YES")) {
+                formatter->setUseApproxTMerc(true);
             }
         }
         obj->lastPROJString = exportable->exportToPROJString(formatter.get());
@@ -1895,6 +1895,8 @@ PROJ_STRING_LIST proj_get_authorities_from_database(PJ_CONTEXT *ctx) {
  *
  * @return a NULL terminated list of NUL-terminated strings that must be
  * freed with proj_string_list_destroy(), or NULL in case of error.
+ *
+ * @see proj_get_crs_info_list_from_database()
  */
 PROJ_STRING_LIST proj_get_codes_from_database(PJ_CONTEXT *ctx,
                                               const char *auth_name,
@@ -1932,6 +1934,215 @@ void proj_string_list_destroy(PROJ_STRING_LIST list) {
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instanciate a default set of parameters to be used by
+ * proj_get_crs_list().
+ *
+ * @return a new object to free with proj_get_crs_list_parameters_destroy() */
+PROJ_CRS_LIST_PARAMETERS *proj_get_crs_list_parameters_create() {
+    auto ret = new (std::nothrow) PROJ_CRS_LIST_PARAMETERS();
+    if (ret) {
+        ret->types = nullptr;
+        ret->typesCount = 0;
+        ret->crs_area_of_use_contains_bbox = TRUE;
+        ret->bbox_valid = FALSE;
+        ret->west_lon_degree = 0.0;
+        ret->south_lat_degree = 0.0;
+        ret->east_lon_degree = 0.0;
+        ret->north_lat_degree = 0.0;
+        ret->allow_deprecated = FALSE;
+    }
+    return ret;
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Destroy an object returned by proj_get_crs_list_parameters_create()
+ */
+void proj_get_crs_list_parameters_destroy(PROJ_CRS_LIST_PARAMETERS *params) {
+    delete params;
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Enumerate CRS objects from the database, taking into account various
+ * criteria.
+ *
+ * The returned object is an array of PROJ_CRS_INFO* pointers, whose last
+ * entry is NULL. This array should be freed with proj_crs_info_list_destroy()
+ *
+ * When no filter parameters are set, this is functionnaly equivalent to
+ * proj_get_crs_info_list_from_database(), instanciating a PJ* object for each
+ * of the proj_create_from_database() and retrieving information with the
+ * various getters. However this function will be much faster.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param auth_name Authority name, used to restrict the search.
+ * Or NULL for all authorities.
+ * @param params Additional criteria, or NULL. If not-NULL, params SHOULD
+ * have been allocated by proj_get_crs_list_parameters_create(), as the
+ * PROJ_CRS_LIST_PARAMETERS structure might grow over time.
+ * @param out_result_count Output parameter pointing to an integer to receive
+ * the size of the result list. Might be NULL
+ * @return an array of PROJ_CRS_INFO* pointers to be freed with
+ * proj_crs_info_list_destroy(), or NULL in case of error.
+ */
+PROJ_CRS_INFO **
+proj_get_crs_info_list_from_database(PJ_CONTEXT *ctx, const char *auth_name,
+                                     const PROJ_CRS_LIST_PARAMETERS *params,
+                                     int *out_result_count) {
+    SANITIZE_CTX(ctx);
+    PROJ_CRS_INFO **ret = nullptr;
+    int i = 0;
+    try {
+        auto factory = AuthorityFactory::create(getDBcontext(ctx),
+                                                auth_name ? auth_name : "");
+        auto list = factory->getCRSInfoList();
+        ret = new PROJ_CRS_INFO *[list.size() + 1];
+        GeographicBoundingBoxPtr bbox;
+        if (params && params->bbox_valid) {
+            bbox = GeographicBoundingBox::create(
+                       params->west_lon_degree, params->south_lat_degree,
+                       params->east_lon_degree, params->north_lat_degree)
+                       .as_nullable();
+        }
+        for (const auto &info : list) {
+            auto type = PJ_TYPE_CRS;
+            if (info.type == AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS) {
+                type = PJ_TYPE_GEOGRAPHIC_2D_CRS;
+            } else if (info.type ==
+                       AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS) {
+                type = PJ_TYPE_GEOGRAPHIC_3D_CRS;
+            } else if (info.type ==
+                       AuthorityFactory::ObjectType::GEOCENTRIC_CRS) {
+                type = PJ_TYPE_GEOCENTRIC_CRS;
+            } else if (info.type ==
+                       AuthorityFactory::ObjectType::PROJECTED_CRS) {
+                type = PJ_TYPE_PROJECTED_CRS;
+            } else if (info.type ==
+                       AuthorityFactory::ObjectType::VERTICAL_CRS) {
+                type = PJ_TYPE_VERTICAL_CRS;
+            } else if (info.type ==
+                       AuthorityFactory::ObjectType::COMPOUND_CRS) {
+                type = PJ_TYPE_COMPOUND_CRS;
+            }
+            if (params && params->typesCount) {
+                bool typeValid = false;
+                for (size_t j = 0; j < params->typesCount; j++) {
+                    if (params->types[j] == type) {
+                        typeValid = true;
+                        break;
+                    } else if (params->types[j] == PJ_TYPE_GEOGRAPHIC_CRS &&
+                               (type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
+                                type == PJ_TYPE_GEOGRAPHIC_3D_CRS)) {
+                        typeValid = true;
+                        break;
+                    } else if (params->types[j] == PJ_TYPE_GEODETIC_CRS &&
+                               (type == PJ_TYPE_GEOCENTRIC_CRS ||
+                                type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
+                                type == PJ_TYPE_GEOGRAPHIC_3D_CRS)) {
+                        typeValid = true;
+                        break;
+                    }
+                }
+                if (!typeValid) {
+                    continue;
+                }
+            }
+            if (params && !params->allow_deprecated && info.deprecated) {
+                continue;
+            }
+            if (params && params->bbox_valid) {
+                if (!info.bbox_valid) {
+                    continue;
+                }
+                if (info.west_lon_degree <= info.east_lon_degree &&
+                    params->west_lon_degree <= params->east_lon_degree) {
+                    if (params->crs_area_of_use_contains_bbox) {
+                        if (params->west_lon_degree < info.west_lon_degree ||
+                            params->east_lon_degree > info.east_lon_degree ||
+                            params->south_lat_degree < info.south_lat_degree ||
+                            params->north_lat_degree > info.north_lat_degree) {
+                            continue;
+                        }
+                    } else {
+                        if (info.east_lon_degree < params->west_lon_degree ||
+                            info.west_lon_degree > params->east_lon_degree ||
+                            info.north_lat_degree < params->south_lat_degree ||
+                            info.south_lat_degree > params->north_lat_degree) {
+                            continue;
+                        }
+                    }
+                } else {
+                    auto crsExtent = GeographicBoundingBox::create(
+                        info.west_lon_degree, info.south_lat_degree,
+                        info.east_lon_degree, info.north_lat_degree);
+                    if (params->crs_area_of_use_contains_bbox) {
+                        if (!crsExtent->contains(NN_NO_CHECK(bbox))) {
+                            continue;
+                        }
+                    } else {
+                        if (!bbox->intersects(crsExtent)) {
+                            continue;
+                        }
+                    }
+                }
+            }
+
+            ret[i] = new PROJ_CRS_INFO;
+            ret[i]->auth_name = pj_strdup(info.authName.c_str());
+            ret[i]->code = pj_strdup(info.code.c_str());
+            ret[i]->name = pj_strdup(info.name.c_str());
+            ret[i]->type = type;
+            ret[i]->deprecated = info.deprecated;
+            ret[i]->bbox_valid = info.bbox_valid;
+            ret[i]->west_lon_degree = info.west_lon_degree;
+            ret[i]->south_lat_degree = info.south_lat_degree;
+            ret[i]->east_lon_degree = info.east_lon_degree;
+            ret[i]->north_lat_degree = info.north_lat_degree;
+            ret[i]->area_name = pj_strdup(info.areaName.c_str());
+            ret[i]->projection_method_name =
+                info.projectionMethodName.empty()
+                    ? nullptr
+                    : pj_strdup(info.projectionMethodName.c_str());
+            i++;
+        }
+        ret[i] = nullptr;
+        if (out_result_count)
+            *out_result_count = i;
+        return ret;
+    } catch (const std::exception &e) {
+        proj_log_error(ctx, __FUNCTION__, e.what());
+        if (ret) {
+            ret[i + 1] = nullptr;
+            proj_crs_info_list_destroy(ret);
+        }
+        if (out_result_count)
+            *out_result_count = 0;
+    }
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Destroy the result returned by
+ * proj_get_crs_info_list_from_database().
+ */
+void proj_crs_info_list_destroy(PROJ_CRS_INFO **list) {
+    if (list) {
+        for (int i = 0; list[i] != nullptr; i++) {
+            pj_dalloc(list[i]->auth_name);
+            pj_dalloc(list[i]->code);
+            pj_dalloc(list[i]->name);
+            pj_dalloc(list[i]->area_name);
+            pj_dalloc(list[i]->projection_method_name);
+            delete list[i];
+        }
+        delete[] list;
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Return the Conversion of a DerivedCRS (such as a ProjectedCRS),
  * or the Transformation from the baseCRS to the hubCRS of a BoundCRS
  *
diff --git a/src/iso19111/common.cpp b/src/iso19111/common.cpp
index 57654d84..bdd836e1 100644
--- a/src/iso19111/common.cpp
+++ b/src/iso19111/common.cpp
@@ -172,21 +172,22 @@ void UnitOfMeasure::_exportToWKT(
 {
     const bool isWKT2 = formatter->version() == WKTFormatter::Version::WKT2;
 
-    if (formatter->forceUNITKeyword() && type() != Type::PARAMETRIC) {
+    const auto l_type = type();
+    if (formatter->forceUNITKeyword() && l_type != Type::PARAMETRIC) {
         formatter->startNode(WKTConstants::UNIT, !codeSpace().empty());
     } else if (!unitType.empty()) {
         formatter->startNode(unitType, !codeSpace().empty());
     } else {
-        if (isWKT2 && type() == Type::LINEAR) {
+        if (isWKT2 && l_type == Type::LINEAR) {
             formatter->startNode(WKTConstants::LENGTHUNIT,
                                  !codeSpace().empty());
-        } else if (isWKT2 && type() == Type::ANGULAR) {
+        } else if (isWKT2 && l_type == Type::ANGULAR) {
             formatter->startNode(WKTConstants::ANGLEUNIT, !codeSpace().empty());
-        } else if (isWKT2 && type() == Type::SCALE) {
+        } else if (isWKT2 && l_type == Type::SCALE) {
             formatter->startNode(WKTConstants::SCALEUNIT, !codeSpace().empty());
-        } else if (isWKT2 && type() == Type::TIME) {
+        } else if (isWKT2 && l_type == Type::TIME) {
             formatter->startNode(WKTConstants::TIMEUNIT, !codeSpace().empty());
-        } else if (isWKT2 && type() == Type::PARAMETRIC) {
+        } else if (isWKT2 && l_type == Type::PARAMETRIC) {
             formatter->startNode(WKTConstants::PARAMETRICUNIT,
                                  !codeSpace().empty());
         } else {
@@ -211,7 +212,7 @@ void UnitOfMeasure::_exportToWKT(
             formatter->addQuotedString(l_name);
         }
         const auto &factor = conversionToSI();
-        if (!isWKT2 || factor != 0.0) {
+        if (!isWKT2 || l_type != Type::TIME || factor != 0.0) {
             // Some TIMEUNIT do not have a conversion factor
             formatter->add(factor);
         }
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 6f9b6283..d36d901f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -2147,11 +2147,19 @@ void ParameterValue::_exportToWKT(io::WKTFormatter *formatter) const {
             // registered linear / angular unit.
             const auto &unitType = unit.type();
             if (unitType == common::UnitOfMeasure::Type::LINEAR) {
-                formatter->add(
-                    l_value.convertToUnit(*(formatter->axisLinearUnit())));
+                const auto &targetUnit = *(formatter->axisLinearUnit());
+                if (targetUnit.conversionToSI() == 0.0) {
+                    throw io::FormattingException(
+                        "cannot convert value to target linear unit");
+                }
+                formatter->add(l_value.convertToUnit(targetUnit));
             } else if (unitType == common::UnitOfMeasure::Type::ANGULAR) {
-                formatter->add(
-                    l_value.convertToUnit(*(formatter->axisAngularUnit())));
+                const auto &targetUnit = *(formatter->axisAngularUnit());
+                if (targetUnit.conversionToSI() == 0.0) {
+                    throw io::FormattingException(
+                        "cannot convert value to target angular unit");
+                }
+                formatter->add(l_value.convertToUnit(targetUnit));
             } else {
                 formatter->add(l_value.getSIValue());
             }
@@ -5373,13 +5381,11 @@ bool Conversion::addWKTExtensionNode(io::WKTFormatter *formatter) const {
         const auto &l_method = method();
         const auto &methodName = l_method->nameStr();
         const int methodEPSGCode = l_method->getEPSGCode();
-        int zone = 0;
-        bool north = true;
-        if (l_method->getPrivate()->projMethodOverride_ == "etmerc" &&
-            !isUTM(zone, north)) {
+        if (l_method->getPrivate()->projMethodOverride_ == "tmerc approx" ||
+            l_method->getPrivate()->projMethodOverride_ == "utm approx") {
             auto projFormatter = io::PROJStringFormatter::create();
             projFormatter->setCRSExport(true);
-            projFormatter->setUseETMercForTMerc(true);
+            projFormatter->setUseApproxTMerc(true);
             formatter->startNode(io::WKTConstants::EXTENSION, false);
             formatter->addQuotedString("PROJ4");
             _exportToPROJString(projFormatter.get());
@@ -5471,17 +5477,21 @@ void Conversion::_exportToPROJString(
     const auto &convName = nameStr();
     bool bConversionDone = false;
     bool bEllipsoidParametersDone = false;
-    bool useETMerc = false;
+    bool useApprox = false;
     if (methodEPSGCode == EPSG_CODE_METHOD_TRANSVERSE_MERCATOR) {
         // Check for UTM
         int zone = 0;
         bool north = true;
-        bool etMercSettingSet = false;
-        useETMerc = formatter->getUseETMercForTMerc(etMercSettingSet) ||
-                    l_method->getPrivate()->projMethodOverride_ == "etmerc";
-        if (isUTM(zone, north) && !(etMercSettingSet && !useETMerc)) {
+        useApprox =
+            formatter->getUseApproxTMerc() ||
+            l_method->getPrivate()->projMethodOverride_ == "tmerc approx" ||
+            l_method->getPrivate()->projMethodOverride_ == "utm approx";
+        if (isUTM(zone, north)) {
             bConversionDone = true;
             formatter->addStep("utm");
+            if( useApprox ) {
+                formatter->addParam("approx");
+            }
             formatter->addParam("zone", zone);
             if (!north) {
                 formatter->addParam("south");
@@ -5654,7 +5664,10 @@ void Conversion::_exportToPROJString(
     if (!bConversionDone) {
         const MethodMapping *mapping = getMapping(l_method.get());
         if (mapping && mapping->proj_name_main) {
-            formatter->addStep(useETMerc ? "etmerc" : mapping->proj_name_main);
+            formatter->addStep(mapping->proj_name_main);
+            if (useApprox) {
+                formatter->addParam("approx");
+            }
             if (mapping->proj_name_aux) {
                 if (internal::starts_with(mapping->proj_name_aux, "axis=")) {
                     bAxisSpecFound = true;
diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp
index 4df1b319..66717f70 100644
--- a/src/iso19111/datum.cpp
+++ b/src/iso19111/datum.cpp
@@ -646,7 +646,8 @@ EllipsoidNNPtr Ellipsoid::createSphere(const util::PropertyMap &properties,
  * @param properties See \ref general_properties.
  * At minimum the name should be defined.
  * @param semiMajorAxisIn the semi-major axis.
- * @param invFlattening the inverse/reverse flattening.
+ * @param invFlattening the inverse/reverse flattening. If set to 0, this will
+ * be considered as a sphere.
  * @param celestialBody Name of the celestial body on which the ellipsoid refers
  * to.
  * @return new Ellipsoid.
@@ -654,8 +655,11 @@ EllipsoidNNPtr Ellipsoid::createSphere(const util::PropertyMap &properties,
 EllipsoidNNPtr Ellipsoid::createFlattenedSphere(
     const util::PropertyMap &properties, const common::Length &semiMajorAxisIn,
     const common::Scale &invFlattening, const std::string &celestialBody) {
-    auto ellipsoid(Ellipsoid::nn_make_shared<Ellipsoid>(
-        semiMajorAxisIn, invFlattening, celestialBody));
+    auto ellipsoid(invFlattening.value() == 0
+                       ? Ellipsoid::nn_make_shared<Ellipsoid>(semiMajorAxisIn,
+                                                              celestialBody)
+                       : Ellipsoid::nn_make_shared<Ellipsoid>(
+                             semiMajorAxisIn, invFlattening, celestialBody));
     ellipsoid->setProperties(properties);
     return ellipsoid;
 }
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 54fc3fd3..81abcdf1 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -77,9 +77,17 @@ namespace io {
 
 //! @cond Doxygen_Suppress
 
-#define GEOG_2D "'geographic 2D'"
-#define GEOG_3D "'geographic 3D'"
-#define GEOCENTRIC "'geocentric'"
+// CRS subtypes
+#define GEOG_2D "geographic 2D"
+#define GEOG_3D "geographic 3D"
+#define GEOCENTRIC "geocentric"
+#define PROJECTED "projected"
+#define VERTICAL "vertical"
+#define COMPOUND "compound"
+
+#define GEOG_2D_SINGLE_QUOTED "'geographic 2D'"
+#define GEOG_3D_SINGLE_QUOTED "'geographic 3D'"
+#define GEOCENTRIC_SINGLE_QUOTED "'geocentric'"
 
 // ---------------------------------------------------------------------------
 
@@ -1051,7 +1059,7 @@ DatabaseContext::getAliasFromOfficialName(const std::string &officialName,
     sql += replaceAll(tableName, "\"", "\"\"");
     sql += "\" WHERE name = ?";
     if (tableName == "geodetic_crs") {
-        sql += " AND type = " GEOG_2D;
+        sql += " AND type = " GEOG_2D_SINGLE_QUOTED;
     }
     auto res = d->run(sql, {officialName});
     if (res.empty()) {
@@ -1485,6 +1493,12 @@ AuthorityFactory::createExtent(const std::string &code) const {
     try {
         const auto &row = res.front();
         const auto &name = row[0];
+        if (row[1].empty()) {
+            auto extent = metadata::Extent::create(
+                util::optional<std::string>(name), {}, {}, {});
+            d->context()->d->cache(code, extent);
+            return extent;
+        }
         double south_lat = c_locale_stod(row[1]);
         double north_lat = c_locale_stod(row[2]);
         double west_lon = c_locale_stod(row[3]);
@@ -1733,7 +1747,7 @@ AuthorityFactory::createEllipsoid(const std::string &code) const {
                 common::Length(c_locale_stod(semi_minor_axis_str), uom), body);
         }
     } catch (const std::exception &ex) {
-        throw buildFactoryException("elllipsoid", code, ex);
+        throw buildFactoryException("ellipsoid", code, ex);
     }
 }
 
@@ -2061,7 +2075,8 @@ AuthorityFactory::createGeodeticCRS(const std::string &code,
                     "deprecated FROM "
                     "geodetic_crs WHERE auth_name = ? AND code = ?");
     if (geographicOnly) {
-        sql += " AND type in (" GEOG_2D "," GEOG_3D ")";
+        sql += " AND type in (" GEOG_2D_SINGLE_QUOTED "," GEOG_3D_SINGLE_QUOTED
+               ")";
     }
     auto res = d->runWithCodeParam(sql, code);
     if (res.empty()) {
@@ -2118,15 +2133,14 @@ AuthorityFactory::createGeodeticCRS(const std::string &code,
 
         auto ellipsoidalCS =
             util::nn_dynamic_pointer_cast<cs::EllipsoidalCS>(cs);
-        if ((type == "geographic 2D" || type == "geographic 3D") &&
-            ellipsoidalCS) {
+        if ((type == GEOG_2D || type == GEOG_3D) && ellipsoidalCS) {
             auto crsRet = crs::GeographicCRS::create(
                 props, datum, NN_NO_CHECK(ellipsoidalCS));
             d->context()->d->cache(cacheKey, crsRet);
             return crsRet;
         }
         auto geocentricCS = util::nn_dynamic_pointer_cast<cs::CartesianCS>(cs);
-        if (type == "geocentric" && geocentricCS) {
+        if (type == GEOCENTRIC && geocentricCS) {
             auto crsRet = crs::GeodeticCRS::create(props, datum,
                                                    NN_NO_CHECK(geocentricCS));
             d->context()->d->cache(cacheKey, crsRet);
@@ -2471,17 +2485,16 @@ AuthorityFactory::createCoordinateReferenceSystem(const std::string &code,
                                            code);
     }
     const auto &type = res.front()[0];
-    if (type == "geographic 2D" || type == "geographic 3D" ||
-        type == "geocentric") {
+    if (type == GEOG_2D || type == GEOG_3D || type == GEOCENTRIC) {
         return createGeodeticCRS(code);
     }
-    if (type == "vertical") {
+    if (type == VERTICAL) {
         return createVerticalCRS(code);
     }
-    if (type == "projected") {
+    if (type == PROJECTED) {
         return createProjectedCRS(code);
     }
-    if (allowCompound && type == "compound") {
+    if (allowCompound && type == COMPOUND) {
         return createCompoundCRS(code);
     }
     throw FactoryException("unhandled CRS type: " + type);
@@ -3776,17 +3789,22 @@ AuthorityFactory::getAuthorityCodes(const ObjectType &type,
         sql = "SELECT code FROM geodetic_crs WHERE ";
         break;
     case ObjectType::GEOCENTRIC_CRS:
-        sql = "SELECT code FROM geodetic_crs WHERE type = " GEOCENTRIC " AND ";
+        sql = "SELECT code FROM geodetic_crs WHERE type "
+              "= " GEOCENTRIC_SINGLE_QUOTED " AND ";
         break;
     case ObjectType::GEOGRAPHIC_CRS:
-        sql = "SELECT code FROM geodetic_crs WHERE type IN (" GEOG_2D
-              "," GEOG_3D ") AND ";
+        sql = "SELECT code FROM geodetic_crs WHERE type IN "
+              "(" GEOG_2D_SINGLE_QUOTED "," GEOG_3D_SINGLE_QUOTED ") AND ";
         break;
     case ObjectType::GEOGRAPHIC_2D_CRS:
-        sql = "SELECT code FROM geodetic_crs WHERE type = " GEOG_2D " AND ";
+        sql =
+            "SELECT code FROM geodetic_crs WHERE type = " GEOG_2D_SINGLE_QUOTED
+            " AND ";
         break;
     case ObjectType::GEOGRAPHIC_3D_CRS:
-        sql = "SELECT code FROM geodetic_crs WHERE type = " GEOG_3D " AND ";
+        sql =
+            "SELECT code FROM geodetic_crs WHERE type = " GEOG_3D_SINGLE_QUOTED
+            " AND ";
         break;
     case ObjectType::VERTICAL_CRS:
         sql = "SELECT code FROM vertical_crs WHERE ";
@@ -3852,6 +3870,107 @@ AuthorityFactory::getDescriptionText(const std::string &code) const {
 
 // ---------------------------------------------------------------------------
 
+/** \brief Return a list of information on CRS objects
+ *
+ * This is functionnaly equivalent of listing the codes from an authority,
+ * instanciating
+ * a CRS object for each of them and getting the information from this CRS
+ * object, but this implementation has much less overhead.
+ *
+ * @throw FactoryException
+ */
+std::list<AuthorityFactory::CRSInfo> AuthorityFactory::getCRSInfoList() const {
+    std::string sql = "SELECT c.auth_name, c.code, c.name, c.type, "
+                      "c.deprecated, "
+                      "a.west_lon, a.south_lat, a.east_lon, a.north_lat, "
+                      "a.name, NULL FROM geodetic_crs c "
+                      "JOIN area a ON "
+                      "c.area_of_use_auth_name = a.auth_name AND "
+                      "c.area_of_use_code = a.code";
+    ListOfParams params;
+    if (d->hasAuthorityRestriction()) {
+        sql += " WHERE c.auth_name = ?";
+        params.emplace_back(d->authority());
+    }
+    sql += " UNION ALL ";
+    sql += "SELECT c.auth_name, c.code, c.name, 'projected', "
+           "c.deprecated, "
+           "a.west_lon, a.south_lat, a.east_lon, a.north_lat, "
+           "a.name, conv.method_name FROM projected_crs c "
+           "JOIN area a ON "
+           "c.area_of_use_auth_name = a.auth_name AND "
+           "c.area_of_use_code = a.code "
+           "LEFT JOIN conversion conv ON "
+           "c.conversion_auth_name = conv.auth_name AND "
+           "c.conversion_code = conv.code";
+    if (d->hasAuthorityRestriction()) {
+        sql += " WHERE c.auth_name = ?";
+        params.emplace_back(d->authority());
+    }
+    sql += " UNION ALL ";
+    sql += "SELECT c.auth_name, c.code, c.name, 'vertical', "
+           "c.deprecated, "
+           "a.west_lon, a.south_lat, a.east_lon, a.north_lat, "
+           "a.name, NULL FROM vertical_crs c "
+           "JOIN area a ON "
+           "c.area_of_use_auth_name = a.auth_name AND "
+           "c.area_of_use_code = a.code";
+    if (d->hasAuthorityRestriction()) {
+        sql += " WHERE c.auth_name = ?";
+        params.emplace_back(d->authority());
+    }
+    sql += " UNION ALL ";
+    sql += "SELECT c.auth_name, c.code, c.name, 'compound', "
+           "c.deprecated, "
+           "a.west_lon, a.south_lat, a.east_lon, a.north_lat, "
+           "a.name, NULL FROM compound_crs c "
+           "JOIN area a ON "
+           "c.area_of_use_auth_name = a.auth_name AND "
+           "c.area_of_use_code = a.code";
+    if (d->hasAuthorityRestriction()) {
+        sql += " WHERE c.auth_name = ?";
+        params.emplace_back(d->authority());
+    }
+    auto sqlRes = d->run(sql, params);
+    std::list<AuthorityFactory::CRSInfo> res;
+    for (const auto &row : sqlRes) {
+        AuthorityFactory::CRSInfo info;
+        info.authName = row[0];
+        info.code = row[1];
+        info.name = row[2];
+        const auto &type = row[3];
+        if (type == GEOG_2D) {
+            info.type = AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS;
+        } else if (type == GEOG_3D) {
+            info.type = AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS;
+        } else if (type == GEOCENTRIC) {
+            info.type = AuthorityFactory::ObjectType::GEOCENTRIC_CRS;
+        } else if (type == PROJECTED) {
+            info.type = AuthorityFactory::ObjectType::PROJECTED_CRS;
+        } else if (type == VERTICAL) {
+            info.type = AuthorityFactory::ObjectType::VERTICAL_CRS;
+        } else if (type == COMPOUND) {
+            info.type = AuthorityFactory::ObjectType::COMPOUND_CRS;
+        }
+        info.deprecated = row[4] == "1";
+        if (row[5].empty()) {
+            info.bbox_valid = false;
+        } else {
+            info.bbox_valid = true;
+            info.west_lon_degree = c_locale_stod(row[5]);
+            info.south_lat_degree = c_locale_stod(row[6]);
+            info.east_lon_degree = c_locale_stod(row[7]);
+            info.north_lat_degree = c_locale_stod(row[8]);
+        }
+        info.areaName = row[9];
+        info.projectionMethodName = row[10];
+        res.emplace_back(info);
+    }
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Gets the official name from a possibly alias name.
  *
  * @param aliasedName Alias name.
@@ -4050,23 +4169,25 @@ AuthorityFactory::createObjectsFromName(
                 break;
             case ObjectType::GEOCENTRIC_CRS:
                 addToListStringWithOR(otherConditions,
-                                      "(table_name = " GEOCENTRIC " AND "
-                                      "type = " GEOCENTRIC ")");
+                                      "(table_name = " GEOCENTRIC_SINGLE_QUOTED
+                                      " AND "
+                                      "type = " GEOCENTRIC_SINGLE_QUOTED ")");
                 break;
             case ObjectType::GEOGRAPHIC_CRS:
                 addToListStringWithOR(otherConditions,
                                       "(table_name = 'geodetic_crs' AND "
-                                      "type IN (" GEOG_2D "," GEOG_3D "))");
+                                      "type IN (" GEOG_2D_SINGLE_QUOTED
+                                      "," GEOG_3D_SINGLE_QUOTED "))");
                 break;
             case ObjectType::GEOGRAPHIC_2D_CRS:
                 addToListStringWithOR(otherConditions,
                                       "(table_name = 'geodetic_crs' AND "
-                                      "type = " GEOG_2D ")");
+                                      "type = " GEOG_2D_SINGLE_QUOTED ")");
                 break;
             case ObjectType::GEOGRAPHIC_3D_CRS:
                 addToListStringWithOR(otherConditions,
                                       "(table_name = 'geodetic_crs' AND "
-                                      "type = " GEOG_3D ")");
+                                      "type = " GEOG_3D_SINGLE_QUOTED ")");
                 break;
             case ObjectType::PROJECTED_CRS:
                 addToListString(tableNameList, "'projected_crs'");
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 3517c225..60c28201 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -4746,8 +4746,7 @@ struct PROJStringFormatter::Private {
     bool omitProjLongLatIfPossible_ = false;
     bool omitZUnitConversion_ = false;
     DatabaseContextPtr dbContext_{};
-    bool useETMercForTMerc_ = false;
-    bool useETMercForTMercSet_ = false;
+    bool useApproxTMerc_ = false;
     bool addNoDefs_ = true;
     bool coordOperationOptimizations_ = false;
     bool crsExport_ = false;
@@ -4803,11 +4802,9 @@ PROJStringFormatter::create(Convention conventionIn,
 
 // ---------------------------------------------------------------------------
 
-/** \brief Set whether Extended Transverse Mercator (etmerc) should be used
- * instead of tmerc */
-void PROJStringFormatter::setUseETMercForTMerc(bool flag) {
-    d->useETMercForTMerc_ = flag;
-    d->useETMercForTMercSet_ = true;
+/** \brief Set whether approximate Transverse Mercator or UTM should be used */
+void PROJStringFormatter::setUseApproxTMerc(bool flag) {
+    d->useApproxTMerc_ = flag;
 }
 
 // ---------------------------------------------------------------------------
@@ -5256,9 +5253,8 @@ PROJStringFormatter::Convention PROJStringFormatter::convention() const {
 
 // ---------------------------------------------------------------------------
 
-bool PROJStringFormatter::getUseETMercForTMerc(bool &settingSetOut) const {
-    settingSetOut = d->useETMercForTMercSet_;
-    return d->useETMercForTMerc_;
+bool PROJStringFormatter::getUseApproxTMerc() const {
+    return d->useApproxTMerc_;
 }
 
 // ---------------------------------------------------------------------------
@@ -5285,6 +5281,9 @@ PROJStringSyntaxParser(const std::string &projString, std::vector<Step> &steps,
     const char *c_str = projString.c_str();
     std::vector<std::string> tokens;
 
+    bool hasProj = false;
+    bool hasInit = false;
+    bool hasPipeline = false;
     {
         size_t i = 0;
         while (true) {
@@ -5310,24 +5309,29 @@ PROJStringSyntaxParser(const std::string &projString, std::vector<Step> &steps,
                 }
                 token += c_str[i];
             }
+            if (in_string) {
+                throw ParsingException("Unbalanced double quote");
+            }
             if (token.empty()) {
                 break;
             }
+            if (!hasPipeline &&
+                (token == "proj=pipeline" || token == "+proj=pipeline")) {
+                hasPipeline = true;
+            } else if (!hasProj && (starts_with(token, "proj=") ||
+                                    starts_with(token, "+proj="))) {
+                hasProj = true;
+            } else if (!hasInit && (starts_with(token, "init=") ||
+                                    starts_with(token, "+init="))) {
+                hasInit = true;
+            }
             tokens.emplace_back(token);
         }
     }
 
     bool prevWasTitle = false;
 
-    if (projString.find("proj=pipeline") == std::string::npos) {
-        const bool hasProj = projString.find("proj=") == 0 ||
-                             projString.find("+proj=") == 0 ||
-                             projString.find(" proj=") != std::string::npos ||
-                             projString.find(" +proj=") != std::string::npos;
-        const bool hasInit = projString.find("init=") == 0 ||
-                             projString.find("+init=") == 0 ||
-                             projString.find(" init=") != std::string::npos ||
-                             projString.find(" +init=") != std::string::npos;
+    if (!hasPipeline) {
         if (hasProj || hasInit) {
             steps.push_back(Step());
         }
@@ -7073,8 +7077,10 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                                 : UnitOfMeasure::NONE)));
         }
 
-        if (step.name == "etmerc") {
-            methodMap.set("proj_method", "etmerc");
+        if (step.name == "tmerc" && hasParamValue(step, "approx")) {
+            methodMap.set("proj_method", "tmerc approx");
+        } else if (step.name == "utm" && hasParamValue(step, "approx")) {
+            methodMap.set("proj_method", "utm approx");
         }
 
         conv = Conversion::create(mapWithUnknownName, methodMap, parameters,
@@ -7447,6 +7453,13 @@ static const metadata::ExtentPtr &getExtent(const crs::CRS *crs) {
  */
 BaseObjectNNPtr
 PROJStringParser::createFromPROJString(const std::string &projString) {
+
+    // In some abnormal situations involving init=epsg:XXXX syntax, we could
+    // have infinite loop
+    if (d->ctx_ && d->ctx_->curStringInCreateFromPROJString == projString) {
+        throw ParsingException("invalid PROJ string");
+    }
+
     d->steps_.clear();
     d->title_.clear();
     d->globalParamValues_.clear();
@@ -7755,11 +7768,17 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
         proj_log_func(pj_context, &logger, Logger::log);
         proj_context_use_proj4_init_rules(pj_context, d->usePROJ4InitRules_);
     }
+    if (d->ctx_) {
+        d->ctx_->curStringInCreateFromPROJString = projString;
+    }
     auto pj = pj_create_internal(
         pj_context, (projString.find("type=crs") != std::string::npos
                          ? projString + " +disable_grid_presence_check"
                          : projString)
                         .c_str());
+    if (d->ctx_) {
+        d->ctx_->curStringInCreateFromPROJString.clear();
+    }
     valid = pj != nullptr;
 
     // Remove parameters not understood by PROJ.
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index b595d87e..f28a1d68 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -117,7 +117,6 @@ SET(SRC_LIBPROJ_PROJECTIONS
         projections/wag7.cpp
         projections/lcca.cpp
         projections/geos.cpp
-        projections/etmerc.cpp
         projections/boggs.cpp
         projections/collg.cpp
         projections/comill.cpp
@@ -201,15 +200,15 @@ SET(SRC_LIBPROJ_ISO19111
 
 SET(SRC_LIBPROJ_CORE
         pj_list.h proj_internal.h proj_math.h
-        aasincos.cpp adjlon.cpp bch2bps.cpp bchgen.cpp
-        biveval.cpp dmstor.cpp mk_cheby.cpp auth.cpp
+        aasincos.cpp adjlon.cpp
+        dmstor.cpp auth.cpp
         deriv.cpp ell_set.cpp ellps.cpp errno.cpp
         factors.cpp fwd.cpp init.cpp inv.cpp
         list.cpp malloc.cpp mlfn.cpp msfn.cpp proj_mdist.cpp
         open_lib.cpp param.cpp phi2.cpp pr_list.cpp
         qsfn.cpp strerrno.cpp
         tsfn.cpp units.cpp ctx.cpp log.cpp zpoly1.cpp rtodms.cpp
-        vector1.cpp release.cpp gauss.cpp
+        release.cpp gauss.cpp
         fileapi.cpp
         gc_reader.cpp gridcatalog.cpp
         nad_cvt.cpp nad_init.cpp nad_intr.cpp
diff --git a/src/mk_cheby.cpp b/src/mk_cheby.cpp
deleted file mode 100644
index 0f3b97ed..00000000
--- a/src/mk_cheby.cpp
+++ /dev/null
@@ -1,193 +0,0 @@
-#include "proj.h"
-#include "proj_internal.h"
-static void /* sum coefficients less than res */
-eval(PJ_UV **w, int nu, int nv, double res, PJ_UV *resid) {
-    int i, j;
-    double ab;
-    PJ_UV *s;
-
-    resid->u = resid->v = 0.;
-    for (i = 0; i < nu; ++i) {
-        s = w[i];
-        for (j = 0; j < nv; ++j) {
-            if ((ab = fabs(s->u)) < res)
-                resid->u += ab;
-            if ((ab = fabs(s->v)) < res)
-                resid->v += ab;
-            ++s;
-        }
-    }
-}
-static Tseries * /* create power series structure */
-makeT(int nru, int nrv) {
-    Tseries *T;
-    int i;
-
-    if (!(T = (Tseries *)pj_malloc(sizeof(Tseries))))
-        return nullptr;
-    if (!(T->cu = (struct PW_COEF *)pj_malloc(sizeof(struct PW_COEF) * nru))) {
-        pj_dalloc(T);
-        return nullptr;
-    }
-    if (!(T->cv = (struct PW_COEF *)pj_malloc(sizeof(struct PW_COEF) * nrv))) {
-        pj_dalloc(T->cu);
-        pj_dalloc(T);
-        return nullptr;
-    }
-
-    for (i = 0; i < nru; ++i)
-        T->cu[i].c = nullptr;
-    for (i = 0; i < nrv; ++i)
-        T->cv[i].c = nullptr;
-    return T;
-}
-Tseries *
-mk_cheby(PJ_UV a, PJ_UV b, double res, PJ_UV *resid, PJ_UV (*func)(PJ_UV), 
-         int nu, int nv, int power) {
-    int j, i, nru, nrv, *ncu, *ncv;
-    Tseries *T = nullptr;
-    PJ_UV **w;
-    double cutres;
-
-    if (!(w = (PJ_UV **)vector2(nu, nv, sizeof(PJ_UV))))
-        return nullptr;
-    if (!(ncu = (int *)vector1(nu + nv, sizeof(int)))) {
-        freev2((void **)w, nu);
-        return nullptr;
-    }
-    ncv = ncu + nu;
-    if (!bchgen(a, b, nu, nv, w, func)) {
-        PJ_UV *s;
-        double ab, *p;
-
-        /* analyse coefficients and adjust until residual OK */
-        cutres = res;
-        for (i = 4; i ; --i) {
-            eval(w, nu, nv, cutres, resid);
-            if (resid->u < res && resid->v < res)
-                break;
-            cutres *= 0.5;
-        }
-        if (i <= 0) /* warn of too many tries */
-            resid->u = - resid->u;
-        /* apply cut resolution and set pointers */
-        nru = nrv = 0;
-        for (j = 0; j < nu; ++j) {
-            ncu[j] = ncv[j] = 0; /* clear column maxes */
-            s = w[j];
-            for (i = 0; i < nv; ++i) {
-                if ((ab = fabs(s->u)) < cutres) /* < resolution ? */
-                    s->u = 0.;		/* clear coefficient */
-                else
-                    ncu[j] = i + 1;	/* update column max */
-                if ((ab = fabs(s->v)) < cutres) /* same for v coef's */
-                    s->v = 0.;
-                else
-                    ncv[j] = i + 1;
-                ++s;
-            }
-            if (ncu[j]) nru = j + 1;	/* update row max */
-            if (ncv[j]) nrv = j + 1;
-        }
-        if (power) { /* convert to bivariate power series */
-            if (!bch2bps(a, b, w, nu, nv))
-                goto error;
-            /* possible change in some row counts, so readjust */
-            nru = nrv = 0;
-            for (j = 0; j < nu; ++j) {
-                ncu[j] = ncv[j] = 0; /* clear column maxes */
-                s = w[j];
-                for (i = 0; i < nv; ++i) {
-                    if (s->u != 0.0)
-                        ncu[j] = i + 1;	/* update column max */
-                    if (s->v != 0.0)
-                        ncv[j] = i + 1;
-                    ++s;
-                }
-                if (ncu[j]) nru = j + 1;	/* update row max */
-                if (ncv[j]) nrv = j + 1;
-            }
-            if ((T = makeT(nru, nrv)) != nullptr ) {
-                T->a = a;
-                T->b = b;
-                T->mu = nru - 1;
-                T->mv = nrv - 1;
-                T->power = 1;
-                for (i = 0; i < nru; ++i) /* store coefficient rows for u */
-                {
-                    if ((T->cu[i].m = ncu[i]) != 0)
-                    {
-                        if ((p = T->cu[i].c =
-                             (double *)pj_malloc(sizeof(double) * ncu[i])))
-                            for (j = 0; j < ncu[i]; ++j)
-                                *p++ = (w[i] + j)->u;
-                        else
-                            goto error;
-                    }
-                }
-                for (i = 0; i < nrv; ++i) /* same for v */
-                {
-                    if ((T->cv[i].m = ncv[i]) != 0)
-                    {
-                        if ((p = T->cv[i].c =
-                             (double *)pj_malloc(sizeof(double) * ncv[i])))
-                            for (j = 0; j < ncv[i]; ++j)
-                                *p++ = (w[i] + j)->v;
-                        else
-                            goto error;
-                    }
-                }
-            }
-        } else if ((T = makeT(nru, nrv)) != nullptr) {
-            /* else make returned Chebyshev coefficient structure */
-            T->mu = nru - 1; /* save row degree */
-            T->mv = nrv - 1;
-            T->a.u = a.u + b.u; /* set argument scaling */
-            T->a.v = a.v + b.v;
-            T->b.u = 1. / (b.u - a.u);
-            T->b.v = 1. / (b.v - a.v);
-            T->power = 0;
-            for (i = 0; i < nru; ++i) /* store coefficient rows for u */
-            {
-                if ((T->cu[i].m = ncu[i]) != 0) 
-                {
-                    if ((p = T->cu[i].c =
-                         (double *)pj_malloc(sizeof(double) * ncu[i])))
-                        for (j = 0; j < ncu[i]; ++j)
-                            *p++ = (w[i] + j)->u;
-                    else
-                        goto error;
-                }
-            }
-            for (i = 0; i < nrv; ++i) /* same for v */
-            {
-                if ((T->cv[i].m = ncv[i]) != 0)
-                {
-                    if ((p = T->cv[i].c =
-                         (double *)pj_malloc(sizeof(double) * ncv[i])))
-                        for (j = 0; j < ncv[i]; ++j)
-                            *p++ = (w[i] + j)->v;
-                    else
-                        goto error;
-                }
-            }
-        } else
-            goto error;
-    }
-    goto gohome;
-  error:
-    if (T) { /* pj_dalloc up possible allocations */
-        for (i = 0; i <= T->mu; ++i)
-            if (T->cu[i].c)
-                pj_dalloc(T->cu[i].c);
-        for (i = 0; i <= T->mv; ++i)
-            if (T->cv[i].c)
-                pj_dalloc(T->cv[i].c);
-        pj_dalloc(T);
-    }
-    T = nullptr;
-  gohome:
-    freev2((void **) w, nu);
-    pj_dalloc(ncu);
-    return T;
-}
diff --git a/src/nad_list.h b/src/nad_list.h
deleted file mode 100644
index f82a2ab7..00000000
--- a/src/nad_list.h
+++ /dev/null
@@ -1,6 +0,0 @@
-/* projection list for program nad2nad */
-PROJ_HEAD(lcc, "Lambert Conformal Conic")
-PROJ_HEAD(omerc, "Oblique Mercator")
-PROJ_HEAD(poly, "Polyconic (American)")
-PROJ_HEAD(tmerc, "Transverse Mercator")
-PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)")
diff --git a/src/proj.h b/src/proj.h
index 552f4b26..f25c3228 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -645,6 +645,74 @@ typedef enum
     PJ_CS_TYPE_TEMPORALMEASURE
 } PJ_COORDINATE_SYSTEM_TYPE;
 
+/** \brief Structure given overall description of a CRS.
+ *
+ * This structure may grow over time, and should not be directly allocated by
+ * client code.
+ */
+typedef struct
+{
+    /** Authority name. */
+    char* auth_name;
+    /** Object code. */
+    char* code;
+    /** Object name. */
+    char* name;
+    /** Object type. */
+    PJ_TYPE type;
+    /** Whether the object is deprecated */
+    int deprecated;
+    /** Whereas the west_lon_degree, south_lat_degree, east_lon_degree and
+     * north_lat_degree fields are valid. */
+    int bbox_valid;
+    /** Western-most longitude of the area of use, in degrees. */
+    double west_lon_degree;
+    /** Southern-most latitude of the area of use, in degrees. */
+    double south_lat_degree;
+    /** Eastern-most longitude of the area of use, in degrees. */
+    double east_lon_degree;
+    /** Northern-most latitude of the area of use, in degrees. */
+    double north_lat_degree;
+    /** Name of the area of use. */
+    char* area_name;
+    /** Name of the projection method for a projected CRS. Might be NULL even
+     *for projected CRS in some cases. */
+    char* projection_method_name;
+} PROJ_CRS_INFO;
+
+/** \brief Structure describing optional parameters for proj_get_crs_list();
+ *
+ * This structure may grow over time, and should not be directly allocated by
+ * client code.
+ */
+typedef struct
+{
+    /** Array of allowed object types. Should be NULL if all types are allowed*/
+    const PJ_TYPE* types;
+    /** Size of types. Should be 0 if all types are allowed*/
+    size_t typesCount;
+
+    /** If TRUE and bbox_valid == TRUE, then only CRS whose area of use
+     * entirely contains the specified bounding box will be returned.
+     * If FALSE and bbox_valid == TRUE, then only CRS whose area of use
+     * intersects the specified bounding box will be returned.
+     */
+    int crs_area_of_use_contains_bbox;
+    /** To set to TRUE so that west_lon_degree, south_lat_degree,
+     * east_lon_degree and north_lat_degree fields are taken into account. */
+    int bbox_valid;
+    /** Western-most longitude of the area of use, in degrees. */
+    double west_lon_degree;
+    /** Southern-most latitude of the area of use, in degrees. */
+    double south_lat_degree;
+    /** Eastern-most longitude of the area of use, in degrees. */
+    double east_lon_degree;
+    /** Northern-most latitude of the area of use, in degrees. */
+    double north_lat_degree;
+
+    /** Whether deprecated objects are allowed. Default to FALSE. */
+    int allow_deprecated;
+} PROJ_CRS_LIST_PARAMETERS;
 
 
 /**@}*/
@@ -774,6 +842,19 @@ PROJ_STRING_LIST PROJ_DLL proj_get_codes_from_database(PJ_CONTEXT *ctx,
                                              PJ_TYPE type,
                                              int allow_deprecated);
 
+PROJ_CRS_LIST_PARAMETERS PROJ_DLL *proj_get_crs_list_parameters_create(void);
+
+void PROJ_DLL proj_get_crs_list_parameters_destroy(
+                                        PROJ_CRS_LIST_PARAMETERS* params);
+
+PROJ_CRS_INFO PROJ_DLL **proj_get_crs_info_list_from_database(
+                                      PJ_CONTEXT *ctx,
+                                      const char *auth_name,
+                                      const PROJ_CRS_LIST_PARAMETERS* params,
+                                      int *out_result_count);
+
+void PROJ_DLL proj_crs_info_list_destroy(PROJ_CRS_INFO** list);
+
 /* ------------------------------------------------------------------------- */
 
 
diff --git a/src/proj_internal.h b/src/proj_internal.h
index f5196939..9ffcc2b3 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -653,7 +653,7 @@ struct FACTORS {
 #define PJD_ERR_LAT_0_OR_ALPHA_EQ_90    -33
 #define PJD_ERR_ELLIPSOID_USE_REQUIRED  -34
 #define PJD_ERR_INVALID_UTM_ZONE        -35
-#define PJD_ERR_TCHEBY_VAL_OUT_OF_RANGE -36
+/* -36 no longer used */
 #define PJD_ERR_FAILED_TO_FIND_PROJ     -37
 #define PJD_ERR_FAILED_TO_LOAD_GRID     -38
 #define PJD_ERR_INVALID_M_OR_N          -39
@@ -677,8 +677,8 @@ struct FACTORS {
 #define PJD_ERR_TOO_MANY_INITS          -57
 #define PJD_ERR_INVALID_ARG             -58
 #define PJD_ERR_INCONSISTENT_UNIT       -59
-/* NOTE: Remember to update pj_strerrno.c and transient_error in */
-/* pj_transform.c when adding new value */
+/* NOTE: Remember to update src/strerrno.cpp, src/apps/gie.cpp and transient_error in */
+/* src/transform.cpp when adding new value */
 
 struct projFileAPI_t;
 
@@ -702,6 +702,8 @@ struct projCtx_t {
     const char* (*file_finder) (PJ_CONTEXT *, const char*, void* user_data) = nullptr;
     void* file_finder_user_data = nullptr;
 
+    std::string curStringInCreateFromPROJString{};
+
     projCtx_t() = default;
     projCtx_t(const projCtx_t&);
     ~projCtx_t();
@@ -852,30 +854,6 @@ COMPLEX pj_zpolyd1(COMPLEX, const COMPLEX *, int, COMPLEX *);
 int pj_deriv(PJ_LP, double, const PJ *, struct DERIVS *);
 int pj_factors(PJ_LP, const PJ *, double, struct FACTORS *);
 
-struct PW_COEF {    /* row coefficient structure */
-    int m;          /* number of c coefficients (=0 for none) */
-    double *c;      /* power coefficients */
-};
-
-/* Approximation structures and procedures */
-typedef struct {    /* Chebyshev or Power series structure */
-    PJ_UV a, b;    /* power series range for evaluation */
-                    /* or Chebyshev argument shift/scaling */
-    struct PW_COEF *cu, *cv;
-    int mu, mv;     /* maximum cu and cv index (+1 for count) */
-    int power;      /* != 0 if power series, else Chebyshev */
-} Tseries;
-
-Tseries PROJ_DLL *mk_cheby(PJ_UV, PJ_UV, double, PJ_UV *, PJ_UV (*)(PJ_UV), int, int, int);
-PJ_UV   bpseval(PJ_UV, Tseries *);
-PJ_UV   bcheval(PJ_UV, Tseries *);
-PJ_UV   biveval(PJ_UV, Tseries *);
-void    *vector1(int, int);
-void   **vector2(int, int, int);
-void     freev2(void **v, int nrows);
-int      bchgen(PJ_UV, PJ_UV, int, int, PJ_UV **, PJ_UV(*)(PJ_UV));
-int      bch2bps(PJ_UV, PJ_UV, PJ_UV **, int, int);
-
 /* nadcon related protos */
 PJ_LP             nad_intr(PJ_LP, struct CTABLE *);
 PJ_LP             nad_cvt(PJ_LP, int, struct CTABLE *);
diff --git a/src/projections/etmerc.cpp b/src/projections/etmerc.cpp
deleted file mode 100644
index e75bc168..00000000
--- a/src/projections/etmerc.cpp
+++ /dev/null
@@ -1,362 +0,0 @@
-/*
-** libproj -- library of cartographic projections
-**
-** Copyright (c) 2008   Gerald I. Evenden
-*/
-
-/*
-** Permission is hereby granted, free of charge, to any person obtaining
-** a copy of this software and associated documentation files (the
-** "Software"), to deal in the Software without restriction, including
-** without limitation the rights to use, copy, modify, merge, publish,
-** distribute, sublicense, and/or sell copies of the Software, and to
-** permit persons to whom the Software is furnished to do so, subject to
-** the following conditions:
-**
-** The above copyright notice and this permission notice shall be
-** included in all copies or substantial portions of the Software.
-**
-** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-
-/* The code in this file is largly based upon procedures:
- *
- * Written by: Knud Poder and Karsten Engsager
- *
- * Based on math from: R.Koenig and K.H. Weise, "Mathematische
- * Grundlagen der hoeheren Geodaesie und Kartographie,
- * Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951.
- *
- * Modified and used here by permission of Reference Networks
- * Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark
- *
-*/
-
-#define PJ_LIB__
-
-#include <errno.h>
-
-#include "proj.h"
-#include "proj_internal.h"
-#include "proj_math.h"
-
-
-namespace { // anonymous namespace
-struct pj_opaque {
-    double    Qn;    /* Merid. quad., scaled to the projection */ \
-    double    Zb;    /* Radius vector in polar coord. systems  */ \
-    double    cgb[6]; /* Constants for Gauss -> Geo lat */ \
-    double    cbg[6]; /* Constants for Geo lat -> Gauss */ \
-    double    utg[6]; /* Constants for transv. merc. -> geo */ \
-    double    gtu[6]; /* Constants for geo -> transv. merc. */
-};
-} // anonymous namespace
-
-PROJ_HEAD(etmerc, "Extended Transverse Mercator")
-    "\n\tCyl, Sph\n\tlat_ts=(0)\nlat_0=(0)";
-PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)")
-    "\n\tCyl, Sph\n\tzone= south";
-
-#define PROJ_ETMERC_ORDER 6
-
-#ifdef _GNU_SOURCE
-    inline
-#endif
-static double gatg(double *p1, int len_p1, double B) {
-    double *p;
-    double h = 0, h1, h2 = 0, cos_2B;
-
-    cos_2B = 2*cos(2*B);
-    p = p1 + len_p1;
-    h1 = *--p;
-    while (p - p1) {
-        h = -h2 + cos_2B*h1 + *--p;
-        h2 = h1;
-        h1 = h;
-    }
-    return (B + h*sin(2*B));
-}
-
-/* Complex Clenshaw summation */
-#ifdef _GNU_SOURCE
-    inline
-#endif
-static double clenS(double *a, int size, double arg_r, double arg_i, double *R, double *I) {
-    double      *p, r, i, hr, hr1, hr2, hi, hi1, hi2;
-    double      sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i;
-
-    /* arguments */
-    p = a + size;
-#ifdef _GNU_SOURCE
-    sincos(arg_r, &sin_arg_r, &cos_arg_r);
-#else
-    sin_arg_r  = sin(arg_r);
-    cos_arg_r  = cos(arg_r);
-#endif
-    sinh_arg_i = sinh(arg_i);
-    cosh_arg_i = cosh(arg_i);
-    r          =  2*cos_arg_r*cosh_arg_i;
-    i          = -2*sin_arg_r*sinh_arg_i;
-
-    /* summation loop */
-    hi1 = hr1 = hi = 0;
-    hr = *--p;
-    for (; a - p;) {
-        hr2 = hr1;
-        hi2 = hi1;
-        hr1 = hr;
-        hi1 = hi;
-        hr  = -hr2 + r*hr1 - i*hi1 + *--p;
-        hi  = -hi2 + i*hr1 + r*hi1;
-    }
-
-    r   = sin_arg_r*cosh_arg_i;
-    i   = cos_arg_r*sinh_arg_i;
-    *R  = r*hr - i*hi;
-    *I  = r*hi + i*hr;
-    return *R;
-}
-
-
-/* Real Clenshaw summation */
-static double clens(double *a, int size, double arg_r) {
-    double      *p, r, hr, hr1, hr2, cos_arg_r;
-
-    p = a + size;
-    cos_arg_r  = cos(arg_r);
-    r          =  2*cos_arg_r;
-
-    /* summation loop */
-    hr1 = 0;
-    hr = *--p;
-    for (; a - p;) {
-        hr2 = hr1;
-        hr1 = hr;
-        hr  = -hr2 + r*hr1 + *--p;
-    }
-    return sin (arg_r)*hr;
-}
-
-
-
-static PJ_XY e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
-    PJ_XY xy = {0.0,0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-    double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
-    double Cn = lp.phi, Ce = lp.lam;
-
-    /* ell. LAT, LNG -> Gaussian LAT, LNG */
-    Cn  = gatg (Q->cbg, PROJ_ETMERC_ORDER, Cn);
-    /* Gaussian LAT, LNG -> compl. sph. LAT */
-#ifdef _GNU_SOURCE
-    sincos (Cn, &sin_Cn, &cos_Cn);
-    sincos (Ce, &sin_Ce, &cos_Ce);
-#else
-    sin_Cn = sin (Cn);
-    cos_Cn = cos (Cn);
-    sin_Ce = sin (Ce);
-    cos_Ce = cos (Ce);
-#endif
-
-    Cn     = atan2 (sin_Cn, cos_Ce*cos_Cn);
-    Ce     = atan2 (sin_Ce*cos_Cn,  hypot (sin_Cn, cos_Cn*cos_Ce));
-
-    /* compl. sph. N, E -> ell. norm. N, E */
-    Ce  = asinh ( tan (Ce) );     /* Replaces: Ce  = log(tan(FORTPI + Ce*0.5)); */
-    Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
-    Ce += dCe;
-    if (fabs (Ce) <= 2.623395162778) {
-        xy.y  = Q->Qn * Cn + Q->Zb;  /* Northing */
-        xy.x  = Q->Qn * Ce;          /* Easting  */
-    } else
-        xy.x = xy.y = HUGE_VAL;
-    return xy;
-}
-
-
-
-static PJ_LP e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
-    PJ_LP lp = {0.0,0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-    double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
-    double Cn = xy.y, Ce = xy.x;
-
-    /* normalize N, E */
-    Cn = (Cn - Q->Zb)/Q->Qn;
-    Ce = Ce/Q->Qn;
-
-    if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */
-        /* norm. N, E -> compl. sph. LAT, LNG */
-        Cn += clenS(Q->utg, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
-        Ce += dCe;
-        Ce = atan (sinh (Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - FORTPI); */
-        /* compl. sph. LAT -> Gaussian LAT, LNG */
-#ifdef _GNU_SOURCE
-        sincos (Cn, &sin_Cn, &cos_Cn);
-        sincos (Ce, &sin_Ce, &cos_Ce);
-#else
-        sin_Cn = sin (Cn);
-        cos_Cn = cos (Cn);
-        sin_Ce = sin (Ce);
-        cos_Ce = cos (Ce);
-#endif
-        Ce     = atan2 (sin_Ce, cos_Ce*cos_Cn);
-        Cn     = atan2 (sin_Cn*cos_Ce,  hypot (sin_Ce, cos_Ce*cos_Cn));
-        /* Gaussian LAT, LNG -> ell. LAT, LNG */
-        lp.phi = gatg (Q->cgb,  PROJ_ETMERC_ORDER, Cn);
-        lp.lam = Ce;
-    }
-    else
-        lp.phi = lp.lam = HUGE_VAL;
-    return lp;
-}
-
-
-static PJ *setup(PJ *P) { /* general initialization */
-    double f, n, np, Z;
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-
-    if (P->es <= 0) {
-        return pj_default_destructor(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
-    }
-
-    /* flattening */
-    f = P->es / (1 + sqrt (1 -  P->es)); /* Replaces: f = 1 - sqrt(1-P->es); */
-
-    /* third flattening */
-    np = n = f/(2 - f);
-
-    /* COEF. OF TRIG SERIES GEO <-> GAUSS */
-    /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */
-    /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */
-    /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */
-
-    Q->cgb[0] = n*( 2 + n*(-2/3.0  + n*(-2      + n*(116/45.0 + n*(26/45.0 +
-                n*(-2854/675.0 ))))));
-    Q->cbg[0] = n*(-2 + n*( 2/3.0  + n*( 4/3.0  + n*(-82/45.0 + n*(32/45.0 +
-                n*( 4642/4725.0))))));
-    np     *= n;
-    Q->cgb[1] = np*(7/3.0 + n*( -8/5.0  + n*(-227/45.0 + n*(2704/315.0 +
-                n*( 2323/945.0)))));
-    Q->cbg[1] = np*(5/3.0 + n*(-16/15.0 + n*( -13/9.0  + n*( 904/315.0 +
-                n*(-1522/945.0)))));
-    np     *= n;
-    /* n^5 coeff corrected from 1262/105 -> -1262/105 */
-    Q->cgb[2] = np*( 56/15.0  + n*(-136/35.0 + n*(-1262/105.0 +
-                n*( 73814/2835.0))));
-    Q->cbg[2] = np*(-26/15.0  + n*(  34/21.0 + n*(    8/5.0   +
-                n*(-12686/2835.0))));
-    np     *= n;
-    /* n^5 coeff corrected from 322/35 -> 332/35 */
-    Q->cgb[3] = np*(4279/630.0 + n*(-332/35.0 + n*(-399572/14175.0)));
-    Q->cbg[3] = np*(1237/630.0 + n*( -12/5.0  + n*( -24832/14175.0)));
-    np     *= n;
-    Q->cgb[4] = np*(4174/315.0 + n*(-144838/6237.0 ));
-    Q->cbg[4] = np*(-734/315.0 + n*( 109598/31185.0));
-    np     *= n;
-    Q->cgb[5] = np*(601676/22275.0 );
-    Q->cbg[5] = np*(444337/155925.0);
-
-    /* Constants of the projections */
-    /* Transverse Mercator (UTM, ITM, etc) */
-    np = n*n;
-    /* Norm. mer. quad, K&W p.50 (96), p.19 (38b), p.5 (2) */
-    Q->Qn = P->k0/(1 + n) * (1 + np*(1/4.0 + np*(1/64.0 + np/256.0)));
-    /* coef of trig series */
-    /* utg := ell. N, E -> sph. N, E,  KW p194 (65) */
-    /* gtu := sph. N, E -> ell. N, E,  KW p196 (69) */
-    Q->utg[0] = n*(-0.5  + n*( 2/3.0 + n*(-37/96.0 + n*( 1/360.0 +
-                n*(  81/512.0 + n*(-96199/604800.0))))));
-    Q->gtu[0] = n*( 0.5  + n*(-2/3.0 + n*(  5/16.0 + n*(41/180.0 +
-                n*(-127/288.0 + n*(  7891/37800.0 ))))));
-    Q->utg[1] = np*(-1/48.0 + n*(-1/15.0 + n*(437/1440.0 + n*(-46/105.0 +
-                n*( 1118711/3870720.0)))));
-    Q->gtu[1] = np*(13/48.0 + n*(-3/5.0  + n*(557/1440.0 + n*(281/630.0 +
-                n*(-1983433/1935360.0)))));
-    np      *= n;
-    Q->utg[2] = np*(-17/480.0 + n*(  37/840.0 + n*(  209/4480.0  +
-                n*( -5569/90720.0 ))));
-    Q->gtu[2] = np*( 61/240.0 + n*(-103/140.0 + n*(15061/26880.0 +
-                n*(167603/181440.0))));
-    np      *= n;
-    Q->utg[3] = np*(-4397/161280.0 + n*(  11/504.0 + n*( 830251/7257600.0)));
-    Q->gtu[3] = np*(49561/161280.0 + n*(-179/168.0 + n*(6601661/7257600.0)));
-    np     *= n;
-    Q->utg[4] = np*(-4583/161280.0 + n*(  108847/3991680.0));
-    Q->gtu[4] = np*(34729/80640.0  + n*(-3418889/1995840.0));
-    np     *= n;
-    Q->utg[5] = np*(-20648693/638668800.0);
-    Q->gtu[5] = np*(212378941/319334400.0);
-
-    /* Gaussian latitude value of the origin latitude */
-    Z = gatg (Q->cbg, PROJ_ETMERC_ORDER, P->phi0);
-
-    /* Origin northing minus true northing at the origin latitude */
-    /* i.e. true northing = N - P->Zb                         */
-    Q->Zb  = - Q->Qn*(Z + clens(Q->gtu, PROJ_ETMERC_ORDER, 2*Z));
-    P->inv = e_inverse;
-    P->fwd = e_forward;
-    return P;
-}
-
-
-
-PJ *PROJECTION(etmerc) {
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
-    if (nullptr==Q)
-        return pj_default_destructor (P, ENOMEM);
-    P->opaque = Q;
-   return setup (P);
-}
-
-
-
-/* utm uses etmerc for the underlying projection */
-
-
-PJ *PROJECTION(utm) {
-    long zone;
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
-    if (nullptr==Q)
-        return pj_default_destructor (P, ENOMEM);
-    P->opaque = Q;
-
-    if (P->es == 0.0) {
-        proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
-        return pj_default_destructor(P, ENOMEM);
-    }
-    if (P->lam0 < -1000.0 || P->lam0 > 1000.0) {
-        return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE);
-    }
-
-    P->y0 = pj_param (P->ctx, P->params, "bsouth").i ? 10000000. : 0.;
-    P->x0 = 500000.;
-    if (pj_param (P->ctx, P->params, "tzone").i) /* zone input ? */
-    {
-        zone = pj_param(P->ctx, P->params, "izone").i;
-        if (zone > 0 && zone <= 60)
-            --zone;
-        else {
-            return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE);
-        }
-    }
-    else /* nearest central meridian input */
-    {
-        zone = lround((floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI)));
-        if (zone < 0)
-            zone = 0;
-        else if (zone >= 60)
-            zone = 59;
-    }
-    P->lam0 = (zone + .5) * M_PI / 30. - M_PI;
-    P->k0 = 0.9996;
-    P->phi0 = 0.;
-
-    return setup (P);
-}
diff --git a/src/projections/labrd.cpp b/src/projections/labrd.cpp
index 330c105f..85ab3ddd 100644
--- a/src/projections/labrd.cpp
+++ b/src/projections/labrd.cpp
@@ -108,6 +108,10 @@ PJ *PROJECTION(labrd) {
         return pj_default_destructor (P, ENOMEM);
     P->opaque = Q;
 
+    if (P->phi0 == 0.) {
+        return pj_default_destructor(P, PJD_ERR_LAT_0_IS_ZERO);
+    }
+
     Az = pj_param(P->ctx, P->params, "razi").f;
     sinp = sin(P->phi0);
     t = 1. - P->es * sinp * sinp;
diff --git a/src/projections/tmerc.cpp b/src/projections/tmerc.cpp
index d1938116..c91c5174 100644
--- a/src/projections/tmerc.cpp
+++ b/src/projections/tmerc.cpp
@@ -1,3 +1,15 @@
+/*
+*                   Transverse Mercator implementations
+*
+* In this file two transverse mercator implementations are found. One of Gerald
+* Evenden/John Snyder origin and one of Knud Poder/Karsten Engsager origin. The
+* former is regarded as "approximate" in the following and the latter is "exact".
+* This word choice has been made to distinguish between the two algorithms, where
+* the Evenden/Snyder implementation is the faster, less accurate implementation
+* and the Poder/Engsager algorithm is a slightly slower, but more accurate
+* implementation.
+*/
+
 #define PJ_LIB__
 
 #include <errno.h>
@@ -5,18 +17,32 @@
 
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 
-PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell";
 
+PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell\n\tapprox";
+PROJ_HEAD(etmerc, "Extended Transverse Mercator") "\n\tCyl, Sph";
+PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") "\n\tCyl, Sph\n\tzone= south approx";
 
 namespace { // anonymous namespace
-struct pj_opaque {
+struct pj_opaque_approx {
     double  esp;
     double  ml0;
     double  *en;
 };
+
+struct pj_opaque_exact {
+    double    Qn;     /* Merid. quad., scaled to the projection */
+    double    Zb;     /* Radius vector in polar coord. systems  */
+    double    cgb[6]; /* Constants for Gauss -> Geo lat */
+    double    cbg[6]; /* Constants for Geo lat -> Gauss */
+    double    utg[6]; /* Constants for transv. merc. -> geo */
+    double    gtu[6]; /* Constants for geo -> transv. merc. */
+};
+
 } // anonymous namespace
 
+/* Constants for "approximate" transverse mercator */
 #define EPS10   1.e-10
 #define FC1 1.
 #define FC2 .5
@@ -27,10 +53,18 @@ struct pj_opaque {
 #define FC7 .02380952380952380952
 #define FC8 .01785714285714285714
 
+/* Constant for "exact" transverse mercator */
+#define PROJ_ETMERC_ORDER 6
 
-static PJ_XY e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
+/*****************************************************************************/
+//
+//                  Approximate Transverse Mercator functions
+//
+/*****************************************************************************/
+
+static PJ_XY approx_e_fwd (PJ_LP lp, PJ *P) {
     PJ_XY xy = {0.0, 0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+    struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(P->opaque);
     double al, als, n, cosphi, sinphi, t;
 
     /*
@@ -70,7 +104,7 @@ static PJ_XY e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
 }
 
 
-static PJ_XY s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
+static PJ_XY approx_s_fwd (PJ_LP lp, PJ *P) {
     PJ_XY xy = {0.0,0.0};
     double b, cosphi;
 
@@ -95,7 +129,7 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
         return xy;
     }
 
-    xy.x = static_cast<struct pj_opaque*>(P->opaque)->ml0 * log ((1. + b) / (1. - b));
+    xy.x = static_cast<struct pj_opaque_approx*>(P->opaque)->ml0 * log ((1. + b) / (1. - b));
     xy.y = cosphi * cos (lp.lam) / sqrt (1. - b * b);
 
     b = fabs ( xy.y );
@@ -110,14 +144,14 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
 
     if (lp.phi < 0.)
         xy.y = -xy.y;
-    xy.y = static_cast<struct pj_opaque*>(P->opaque)->esp * (xy.y - P->phi0);
+    xy.y = static_cast<struct pj_opaque_approx*>(P->opaque)->esp * (xy.y - P->phi0);
     return xy;
 }
 
 
-static PJ_LP e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+static PJ_LP approx_e_inv (PJ_XY xy, PJ *P) {
     PJ_LP lp = {0.0,0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+    struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(P->opaque);
     double n, con, cosphi, d, ds, sinphi, t;
 
     lp.phi = pj_inv_mlfn(P->ctx, Q->ml0 + xy.y / P->k0, P->es, Q->en);
@@ -149,13 +183,13 @@ static PJ_LP e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
 }
 
 
-static PJ_LP s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
+static PJ_LP approx_s_inv (PJ_XY xy, PJ *P) {
     PJ_LP lp = {0.0, 0.0};
     double h, g;
 
-    h = exp(xy.x / static_cast<struct pj_opaque*>(P->opaque)->esp);
+    h = exp(xy.x / static_cast<struct pj_opaque_approx*>(P->opaque)->esp);
     g = .5 * (h - 1. / h);
-    h = cos (P->phi0 + xy.y / static_cast<struct pj_opaque*>(P->opaque)->esp);
+    h = cos (P->phi0 + xy.y / static_cast<struct pj_opaque_approx*>(P->opaque)->esp);
     lp.phi = asin(sqrt((1. - h * h) / (1. + g * g)));
 
     /* Make sure that phi is on the correct hemisphere when false northing is used */
@@ -166,45 +200,386 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
 }
 
 
-static PJ *destructor(PJ *P, int errlev) {                       /* Destructor */
+static PJ *destructor_approx(PJ *P, int errlev) {
     if (nullptr==P)
         return nullptr;
 
     if (nullptr==P->opaque)
         return pj_default_destructor(P, errlev);
 
-    pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->en);
+    pj_dealloc (static_cast<struct pj_opaque_approx*>(P->opaque)->en);
     return pj_default_destructor(P, errlev);
 }
 
 
-static PJ *setup(PJ *P) {                   /* general initialization */
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+static PJ *setup_approx(PJ *P) {
+    struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(P->opaque);
+
+    P->destructor = destructor_approx;
+
     if (P->es != 0.0) {
         if (!(Q->en = pj_enfn(P->es)))
             return pj_default_destructor(P, ENOMEM);
 
         Q->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en);
         Q->esp = P->es / (1. - P->es);
-        P->inv = e_inverse;
-        P->fwd = e_forward;
+        P->inv = approx_e_inv;
+        P->fwd = approx_e_fwd;
     } else {
         Q->esp = P->k0;
         Q->ml0 = .5 * Q->esp;
-        P->inv = s_inverse;
-        P->fwd = s_forward;
+        P->inv = approx_s_inv;
+        P->fwd = approx_s_fwd;
     }
     return P;
 }
 
 
+
+/*****************************************************************************/
+//
+//                  Exact Transverse Mercator functions
+//
+//
+// The code in this file is largly based upon procedures:
+//
+// Written by: Knud Poder and Karsten Engsager
+//
+// Based on math from: R.Koenig and K.H. Weise, "Mathematische
+// Grundlagen der hoeheren Geodaesie und Kartographie,
+// Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951.
+//
+// Modified and used here by permission of Reference Networks
+// Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark
+//
+/*****************************************************************************/
+
+/* Helper functios for "exact" transverse mercator */
+#ifdef _GNU_SOURCE
+    inline
+#endif
+static double gatg(double *p1, int len_p1, double B) {
+    double *p;
+    double h = 0, h1, h2 = 0, cos_2B;
+
+    cos_2B = 2*cos(2*B);
+    p = p1 + len_p1;
+    h1 = *--p;
+    while (p - p1) {
+        h = -h2 + cos_2B*h1 + *--p;
+        h2 = h1;
+        h1 = h;
+    }
+    return (B + h*sin(2*B));
+}
+
+/* Complex Clenshaw summation */
+#ifdef _GNU_SOURCE
+    inline
+#endif
+static double clenS(double *a, int size, double arg_r, double arg_i, double *R, double *I) {
+    double      *p, r, i, hr, hr1, hr2, hi, hi1, hi2;
+    double      sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i;
+
+    /* arguments */
+    p = a + size;
+#ifdef _GNU_SOURCE
+    sincos(arg_r, &sin_arg_r, &cos_arg_r);
+#else
+    sin_arg_r  = sin(arg_r);
+    cos_arg_r  = cos(arg_r);
+#endif
+    sinh_arg_i = sinh(arg_i);
+    cosh_arg_i = cosh(arg_i);
+    r          =  2*cos_arg_r*cosh_arg_i;
+    i          = -2*sin_arg_r*sinh_arg_i;
+
+    /* summation loop */
+    hi1 = hr1 = hi = 0;
+    hr = *--p;
+    for (; a - p;) {
+        hr2 = hr1;
+        hi2 = hi1;
+        hr1 = hr;
+        hi1 = hi;
+        hr  = -hr2 + r*hr1 - i*hi1 + *--p;
+        hi  = -hi2 + i*hr1 + r*hi1;
+    }
+
+    r   = sin_arg_r*cosh_arg_i;
+    i   = cos_arg_r*sinh_arg_i;
+    *R  = r*hr - i*hi;
+    *I  = r*hi + i*hr;
+    return *R;
+}
+
+
+/* Real Clenshaw summation */
+static double clens(double *a, int size, double arg_r) {
+    double      *p, r, hr, hr1, hr2, cos_arg_r;
+
+    p = a + size;
+    cos_arg_r  = cos(arg_r);
+    r          =  2*cos_arg_r;
+
+    /* summation loop */
+    hr1 = 0;
+    hr = *--p;
+    for (; a - p;) {
+        hr2 = hr1;
+        hr1 = hr;
+        hr  = -hr2 + r*hr1 + *--p;
+    }
+    return sin (arg_r)*hr;
+}
+
+/* Ellipsoidal, forward */
+static PJ_XY exact_e_fwd (PJ_LP lp, PJ *P) {
+    PJ_XY xy = {0.0,0.0};
+    struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque);
+    double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
+    double Cn = lp.phi, Ce = lp.lam;
+
+    /* ell. LAT, LNG -> Gaussian LAT, LNG */
+    Cn  = gatg (Q->cbg, PROJ_ETMERC_ORDER, Cn);
+    /* Gaussian LAT, LNG -> compl. sph. LAT */
+#ifdef _GNU_SOURCE
+    sincos (Cn, &sin_Cn, &cos_Cn);
+    sincos (Ce, &sin_Ce, &cos_Ce);
+#else
+    sin_Cn = sin (Cn);
+    cos_Cn = cos (Cn);
+    sin_Ce = sin (Ce);
+    cos_Ce = cos (Ce);
+#endif
+
+    Cn     = atan2 (sin_Cn, cos_Ce*cos_Cn);
+    Ce     = atan2 (sin_Ce*cos_Cn,  hypot (sin_Cn, cos_Cn*cos_Ce));
+
+    /* compl. sph. N, E -> ell. norm. N, E */
+    Ce  = asinh ( tan (Ce) );     /* Replaces: Ce  = log(tan(FORTPI + Ce*0.5)); */
+    Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
+    Ce += dCe;
+    if (fabs (Ce) <= 2.623395162778) {
+        xy.y  = Q->Qn * Cn + Q->Zb;  /* Northing */
+        xy.x  = Q->Qn * Ce;          /* Easting  */
+    } else
+        xy.x = xy.y = HUGE_VAL;
+    return xy;
+}
+
+
+/* Ellipsoidal, inverse */
+static PJ_LP exact_e_inv (PJ_XY xy, PJ *P) {
+    PJ_LP lp = {0.0,0.0};
+    struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque);
+    double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
+    double Cn = xy.y, Ce = xy.x;
+
+    /* normalize N, E */
+    Cn = (Cn - Q->Zb)/Q->Qn;
+    Ce = Ce/Q->Qn;
+
+    if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */
+        /* norm. N, E -> compl. sph. LAT, LNG */
+        Cn += clenS(Q->utg, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
+        Ce += dCe;
+        Ce = atan (sinh (Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - FORTPI); */
+        /* compl. sph. LAT -> Gaussian LAT, LNG */
+#ifdef _GNU_SOURCE
+        sincos (Cn, &sin_Cn, &cos_Cn);
+        sincos (Ce, &sin_Ce, &cos_Ce);
+#else
+        sin_Cn = sin (Cn);
+        cos_Cn = cos (Cn);
+        sin_Ce = sin (Ce);
+        cos_Ce = cos (Ce);
+#endif
+        Ce     = atan2 (sin_Ce, cos_Ce*cos_Cn);
+        Cn     = atan2 (sin_Cn*cos_Ce,  hypot (sin_Ce, cos_Ce*cos_Cn));
+        /* Gaussian LAT, LNG -> ell. LAT, LNG */
+        lp.phi = gatg (Q->cgb,  PROJ_ETMERC_ORDER, Cn);
+        lp.lam = Ce;
+    }
+    else
+        lp.phi = lp.lam = HUGE_VAL;
+    return lp;
+}
+
+static PJ *setup_exact(PJ *P) {
+    double f, n, np, Z;
+    struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque);
+
+    if (P->es <= 0) {
+        return pj_default_destructor(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
+    }
+
+    /* flattening */
+    f = P->es / (1 + sqrt (1 -  P->es)); /* Replaces: f = 1 - sqrt(1-P->es); */
+
+    /* third flattening */
+    np = n = f/(2 - f);
+
+    /* COEF. OF TRIG SERIES GEO <-> GAUSS */
+    /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */
+    /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */
+    /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */
+
+    Q->cgb[0] = n*( 2 + n*(-2/3.0  + n*(-2      + n*(116/45.0 + n*(26/45.0 +
+                n*(-2854/675.0 ))))));
+    Q->cbg[0] = n*(-2 + n*( 2/3.0  + n*( 4/3.0  + n*(-82/45.0 + n*(32/45.0 +
+                n*( 4642/4725.0))))));
+    np     *= n;
+    Q->cgb[1] = np*(7/3.0 + n*( -8/5.0  + n*(-227/45.0 + n*(2704/315.0 +
+                n*( 2323/945.0)))));
+    Q->cbg[1] = np*(5/3.0 + n*(-16/15.0 + n*( -13/9.0  + n*( 904/315.0 +
+                n*(-1522/945.0)))));
+    np     *= n;
+    /* n^5 coeff corrected from 1262/105 -> -1262/105 */
+    Q->cgb[2] = np*( 56/15.0  + n*(-136/35.0 + n*(-1262/105.0 +
+                n*( 73814/2835.0))));
+    Q->cbg[2] = np*(-26/15.0  + n*(  34/21.0 + n*(    8/5.0   +
+                n*(-12686/2835.0))));
+    np     *= n;
+    /* n^5 coeff corrected from 322/35 -> 332/35 */
+    Q->cgb[3] = np*(4279/630.0 + n*(-332/35.0 + n*(-399572/14175.0)));
+    Q->cbg[3] = np*(1237/630.0 + n*( -12/5.0  + n*( -24832/14175.0)));
+    np     *= n;
+    Q->cgb[4] = np*(4174/315.0 + n*(-144838/6237.0 ));
+    Q->cbg[4] = np*(-734/315.0 + n*( 109598/31185.0));
+    np     *= n;
+    Q->cgb[5] = np*(601676/22275.0 );
+    Q->cbg[5] = np*(444337/155925.0);
+
+    /* Constants of the projections */
+    /* Transverse Mercator (UTM, ITM, etc) */
+    np = n*n;
+    /* Norm. mer. quad, K&W p.50 (96), p.19 (38b), p.5 (2) */
+    Q->Qn = P->k0/(1 + n) * (1 + np*(1/4.0 + np*(1/64.0 + np/256.0)));
+    /* coef of trig series */
+    /* utg := ell. N, E -> sph. N, E,  KW p194 (65) */
+    /* gtu := sph. N, E -> ell. N, E,  KW p196 (69) */
+    Q->utg[0] = n*(-0.5  + n*( 2/3.0 + n*(-37/96.0 + n*( 1/360.0 +
+                n*(  81/512.0 + n*(-96199/604800.0))))));
+    Q->gtu[0] = n*( 0.5  + n*(-2/3.0 + n*(  5/16.0 + n*(41/180.0 +
+                n*(-127/288.0 + n*(  7891/37800.0 ))))));
+    Q->utg[1] = np*(-1/48.0 + n*(-1/15.0 + n*(437/1440.0 + n*(-46/105.0 +
+                n*( 1118711/3870720.0)))));
+    Q->gtu[1] = np*(13/48.0 + n*(-3/5.0  + n*(557/1440.0 + n*(281/630.0 +
+                n*(-1983433/1935360.0)))));
+    np      *= n;
+    Q->utg[2] = np*(-17/480.0 + n*(  37/840.0 + n*(  209/4480.0  +
+                n*( -5569/90720.0 ))));
+    Q->gtu[2] = np*( 61/240.0 + n*(-103/140.0 + n*(15061/26880.0 +
+                n*(167603/181440.0))));
+    np      *= n;
+    Q->utg[3] = np*(-4397/161280.0 + n*(  11/504.0 + n*( 830251/7257600.0)));
+    Q->gtu[3] = np*(49561/161280.0 + n*(-179/168.0 + n*(6601661/7257600.0)));
+    np     *= n;
+    Q->utg[4] = np*(-4583/161280.0 + n*(  108847/3991680.0));
+    Q->gtu[4] = np*(34729/80640.0  + n*(-3418889/1995840.0));
+    np     *= n;
+    Q->utg[5] = np*(-20648693/638668800.0);
+    Q->gtu[5] = np*(212378941/319334400.0);
+
+    /* Gaussian latitude value of the origin latitude */
+    Z = gatg (Q->cbg, PROJ_ETMERC_ORDER, P->phi0);
+
+    /* Origin northing minus true northing at the origin latitude */
+    /* i.e. true northing = N - P->Zb                         */
+    Q->Zb  = - Q->Qn*(Z + clens(Q->gtu, PROJ_ETMERC_ORDER, 2*Z));
+    P->inv = exact_e_inv;
+    P->fwd = exact_e_fwd;
+    return P;
+}
+
+
+
+
+/*****************************************************************************/
+//
+//                                Operation Setups
+//
+/*****************************************************************************/
+
 PJ *PROJECTION(tmerc) {
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
+    /* exact transverse mercator only exists in ellipsoidal form, */
+    /* use approximate version if +a sphere is requested          */
+    if (pj_param (P->ctx, P->params, "bapprox").i || P->es <= 0) {
+        struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(pj_calloc (1, sizeof (struct pj_opaque_approx)));
+        if (nullptr==Q)
+            return pj_default_destructor (P, ENOMEM);
+
+        P->opaque = Q;
+
+        return setup_approx(P);
+    } else {
+        struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(pj_calloc (1, sizeof (struct pj_opaque_exact)));
+        if (nullptr==Q)
+            return pj_default_destructor (P, ENOMEM);
+        P->opaque = Q;
+       return setup_exact (P);
+    }
+}
+
+
+PJ *PROJECTION(etmerc) {
+    struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(pj_calloc (1, sizeof (struct pj_opaque_exact)));
     if (nullptr==Q)
         return pj_default_destructor (P, ENOMEM);
-
     P->opaque = Q;
-    P->destructor = destructor;
+   return setup_exact (P);
+}
+
 
-    return setup(P);
+/* UTM uses the Poder/Engsager implementation for the underlying projection      */
+/* UNLESS +approx is set in which case the Evenden/Snyder implemenation is used. */
+PJ *PROJECTION(utm) {
+    long zone;
+    if (P->es == 0.0) {
+        proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
+        return pj_default_destructor(P, ENOMEM);
+    }
+    if (P->lam0 < -1000.0 || P->lam0 > 1000.0) {
+        return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE);
+    }
+
+    P->y0 = pj_param (P->ctx, P->params, "bsouth").i ? 10000000. : 0.;
+    P->x0 = 500000.;
+    if (pj_param (P->ctx, P->params, "tzone").i) /* zone input ? */
+    {
+        zone = pj_param(P->ctx, P->params, "izone").i;
+        if (zone > 0 && zone <= 60)
+            --zone;
+        else {
+            return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE);
+        }
+    }
+    else /* nearest central meridian input */
+    {
+        zone = lround((floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI)));
+        if (zone < 0)
+            zone = 0;
+        else if (zone >= 60)
+            zone = 59;
+    }
+    P->lam0 = (zone + .5) * M_PI / 30. - M_PI;
+    P->k0 = 0.9996;
+    P->phi0 = 0.;
+
+    if (pj_param(P->ctx, P->params, "bapprox").i) {
+        struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(pj_calloc (1, sizeof (struct pj_opaque_approx)));
+        if (nullptr==Q)
+            return pj_default_destructor (P, ENOMEM);
+        P->opaque = Q;
+
+        return setup_approx(P);
+    } else {
+        struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(pj_calloc (1, sizeof (struct pj_opaque_exact)));
+        if (nullptr==Q)
+            return pj_default_destructor (P, ENOMEM);
+        P->opaque = Q;
+
+        return setup_exact(P);
+    }
 }
diff --git a/src/strerrno.cpp b/src/strerrno.cpp
index 9f690041..13e9c757 100644
--- a/src/strerrno.cpp
+++ b/src/strerrno.cpp
@@ -44,7 +44,7 @@ pj_err_list[] = {
     "lat_1=lat_2 or lat_1=0 or lat_2=90",                              /* -33 */
     "elliptical usage required",                                       /* -34 */
     "invalid UTM zone number",                                         /* -35 */
-    "arg(s) out of range for Tcheby eval",                             /* -36 */
+    "", /* no longer used */                                           /* -36 */
     "failed to find projection to be rotated",                         /* -37 */
     "failed to load datum shift file",                                 /* -38 */
     "both n & m must be spec'd and > 0",                               /* -39 */
diff --git a/src/vector1.cpp b/src/vector1.cpp
deleted file mode 100644
index fc69f5c3..00000000
--- a/src/vector1.cpp
+++ /dev/null
@@ -1,30 +0,0 @@
-/* make storage for one and two dimensional matricies */
-#include <stdlib.h>
-#include "proj.h"
-#include "proj_internal.h"
-	void * /* one dimension array */
-vector1(int nvals, int size) { return((void *)pj_malloc(size * nvals)); }
-	void /* free 2D array */
-freev2(void **v, int nrows) {
-	if (v) {
-		for (v += nrows; nrows > 0; --nrows)
-			pj_dalloc(*--v);
-		pj_dalloc(v);
-	}
-}
-	void ** /* two dimension array */
-vector2(int nrows, int ncols, int size) {
-	void **s;
-
-	if ((s = (void **)pj_malloc(sizeof(void *) * nrows)) != nullptr) {
-		int rsize, i;
-
-		rsize = size * ncols;
-		for (i = 0; i < nrows; ++i)
-			if (!(s[i] = pj_malloc(rsize))) {
-				freev2(s, i);
-				return (void **)nullptr;
-			}
-	}
-	return s;
-}
diff --git a/test/Makefile.am b/test/Makefile.am
index b65e9785..fa2c1e80 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -7,4 +7,7 @@ EXTRA_DIST = CMakeLists.txt
 # libtestmain_la_LDFLAGS = -no-undefined
 # libtestmain_la_SOURCES = test_main.cpp
 
-SUBDIRS = cli googletest . gie gigs unit
+if !USE_EXTERNAL_GTEST
+MAYBE_GTEST = googletest
+endif
+SUBDIRS = cli $(MAYBE_GTEST) . gie gigs unit
diff --git a/test/cli/Makefile.am b/test/cli/Makefile.am
index 50654968..c0bc0871 100644
--- a/test/cli/Makefile.am
+++ b/test/cli/Makefile.am
@@ -5,6 +5,7 @@ EXEPATH = ../../src
 PROJEXE = $(EXEPATH)/proj
 CS2CSEXE = $(EXEPATH)/cs2cs
 PROJINFOEXE = $(EXEPATH)/projinfo
+CCTEXE = $(EXEPATH)/cct
 
 # PROJ.4 test scripts
 TEST27 = $(THIS_DIR)/test27
@@ -15,12 +16,14 @@ TESTFLAKY = $(THIS_DIR)/testflaky
 TESTDATUMFILE = $(THIS_DIR)/testdatumfile
 TESTIGN = $(THIS_DIR)/testIGNF
 TESTPROJINFO = $(THIS_DIR)/testprojinfo
+TESTCCT = $(THIS_DIR)/testcct
 
 EXTRA_DIST = pj_out27.dist pj_out83.dist td_out.dist \
 		test27 test83 tv_out.dist tf_out.dist \
 		testflaky testvarious testdatumfile testntv2 ntv2_out.dist \
 		testIGNF proj_outIGNF.dist \
 		testprojinfo testprojinfo_out.dist \
+		testcct testcct_out.dist \
 		CMakeLists.txt
 
 testprojinfo-check:
@@ -50,4 +53,7 @@ testntv2-check:
 	  PROJ_LIB=$(DATAPATH) $(TESTNTV2) $(CS2CSEXE) ; \
 	fi
 
-check-local: testprojinfo-check test27-check test83-check testvarious-check testdatumfile-check testign-check testntv2-check
+testcct-check:
+	PROJ_LIB=$(DATAPATH) $(TESTCCT) $(CCTEXE)
+
+check-local: testprojinfo-check test27-check test83-check testvarious-check testdatumfile-check testign-check testntv2-check testcct-check
diff --git a/test/cli/testcct b/test/cli/testcct
new file mode 100755
index 00000000..93749052
--- /dev/null
+++ b/test/cli/testcct
@@ -0,0 +1,51 @@
+:
+# Test cct
+
+TEST_CLI_DIR=`dirname $0`
+DATA_DIR=`dirname $0`/../../data
+EXE=$1
+
+usage()
+{
+    echo "Usage: ${0} <path to 'cct' program>" 
+    echo
+    exit 1
+}
+
+if test -z "${EXE}"; then
+    EXE=../../src/cct
+fi
+
+if test ! -x ${EXE}; then
+    echo "*** ERROR: Can not find '${EXE}' program!"
+    exit 1
+fi
+
+echo "============================================"
+echo "Running ${0} using ${EXE}:"
+echo "============================================"
+
+OUT=testcct_out
+
+rm -f ${OUT}
+
+echo "Testing cct -d 8 +proj=merc +R=1" >> ${OUT}
+echo "90 45" 0 | $EXE -d 8 +proj=merc +R=1 >>${OUT}
+echo "" >>${OUT}
+
+# do 'diff' with distribution results
+echo "diff ${OUT} with testcct_out.dist"
+diff -u ${OUT} ${TEST_CLI_DIR}/testcct_out.dist
+if [ $? -ne 0 ] ; then
+	echo  ""
+	echo "PROBLEMS HAVE OCCURRED"
+	echo "test file ${OUT} saved"
+    echo
+	exit 100
+else
+	echo "TEST OK"
+	echo "test file ${OUT} removed"
+    echo
+	/bin/rm -f ${OUT}
+	exit 0
+fi
diff --git a/test/cli/testcct_out.dist b/test/cli/testcct_out.dist
new file mode 100644
index 00000000..44dd6964
--- /dev/null
+++ b/test/cli/testcct_out.dist
@@ -0,0 +1,3 @@
+Testing cct -d 8 +proj=merc +R=1
+   1.57079633     0.88137359    0.00000000           inf
+
diff --git a/test/cli/testvarious b/test/cli/testvarious
index 7ec50bb3..c1fa61df 100755
--- a/test/cli/testvarious
+++ b/test/cli/testvarious
@@ -239,7 +239,7 @@ EOF
 echo "##############################################################" >> ${OUT}
 echo "Test transverse mercator (#97)" >> ${OUT}
 #
-$EXE +proj=tmerc +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \
+$EXE +proj=tmerc +approx +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \
  +to +proj=latlong +datum=WGS84 \
  -E >>${OUT} <<EOF
 10000 20000
@@ -253,7 +253,7 @@ echo "##############################################################" >> ${OUT}
 echo "Test transverse mercator inverse (#97)" >> ${OUT}
 #
 $EXE +proj=latlong +datum=WGS84 \
- +to +proj=tmerc +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \
+ +to +proj=tmerc +approx +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \
  -E >>${OUT} <<EOF
 0dN 0.000
 15d22'16.108"W	17d52'53.478"N 0.000
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index 148d413d..72e95634 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -352,7 +352,7 @@ Test inverse handling
 10 20	-1384841.19	7581707.88 0.00
 ##############################################################
 Test MGI datum gives expected results (#207)
-16.33 48.20	595710.3732102	5357598.4645755 -44.4951085
+16.33 48.20	595710.3731286	5357598.4645652 -44.4951085
 ##############################################################
 Test omerc sensitivity with locations 90d from origin(#114)
 56.958381652832 72.8798	-9985.16336453	-227.67701050 0.00000000
diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie
index f66b26ec..215971a0 100644
--- a/test/gie/4D-API_cs2cs-style.gie
+++ b/test/gie/4D-API_cs2cs-style.gie
@@ -348,6 +348,14 @@ accept      0 0 1000
 expect      0 0 1
 roundtrip   1
 
+operation   +proj=longlat +a=1 +b=1 +vto_meter=2000/2
+accept      0 0 1000
+expect      0 0 1
+roundtrip   1
+
+operation  +proj=longlat +a=1 +b=1 +vto_meter=1/0
+expect     failure   errno unit_factor_less_than_0
+
 operation   +proj=longlat +a=1 +b=1 +vto_meter=1000 +geoc
 accept      0 0 1000
 expect      0 0 1
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index bc9374a5..930ca71f 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -2243,6 +2243,10 @@ expect  0.498202281 2.000904357
 accept  -200 -100
 expect  0.498202283 1.999095641
 
+-------------------------------------------------------------------------------
+operation +proj=labrd   +ellps=GRS80 +lat_0=0
+accept  0 0
+expect  failure errno lat_0_is_zero
 
 ===============================================================================
 Lambert Azimuthal Equal Area
@@ -5214,6 +5218,16 @@ accept  -200 -100
 expect  -7.490535682 -0.000901935
 
 
+operation +proj=utm +zone=32
+tolerance 0.001 mm
+accept	12 56 0 2000
+expect  687071.43910944  6210141.32674801    0.00000000     2000.0000
+
+operation +proj=utm +zone=32 +approx
+tolerance 0.001 mm
+accept	12 56 0 2000
+expect  687071.43911000  6210141.32675053    0.00000000     2000.0000
+
 ===============================================================================
 van der Grinten (I)
 	Misc Sph
diff --git a/test/gie/ellipsoid.gie b/test/gie/ellipsoid.gie
index 929eb799..ea7f3888 100644
--- a/test/gie/ellipsoid.gie
+++ b/test/gie/ellipsoid.gie
@@ -64,6 +64,9 @@ expect    0 0
 operation proj=merc +es=-1
 expect    failure errno major_axis_not_given
 
+operation proj=merc +R=0
+expect    failure errno major_axis_not_given
+
 operation
 expect    failure
 operation cobra
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 44d01385..d1c978a3 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -496,6 +496,14 @@ expect     failure   errno unit_factor_less_than_0
 operation  proj=utm  ellps=GRS80  zone=32  to_meter=10
 accept     12   55
 expect        69187.5632    609890.7825
+
+operation  proj=utm  ellps=GRS80  zone=32  to_meter=1/0
+expect     failure   errno unit_factor_less_than_0
+
+operation  proj=utm  ellps=GRS80  zone=32  to_meter=2.0/0.2
+accept     12   55
+expect        69187.5632    609890.7825
+
 -------------------------------------------------------------------------------
 
 -------------------------------------------------------------------------------
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index a157f630..86e3484a 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -1,5 +1,16 @@
 # CMake configuration for PROJ unit tests
 
+option(USE_EXTERNAL_GTEST "Compile against external GTest" OFF)
+
+if (USE_EXTERNAL_GTEST)
+
+message(STATUS "Using external GTest")
+find_package(GTest 1.8.0)
+
+else (USE_EXTERNAL_GTEST)
+
+message(STATUS "Using internal GTest")
+
 # FIXME: Deal with our old-school CMakeLists.txt behaving badly
 set(_save_c_flags "${CMAKE_C_FLAGS}")
 set(_save_cxx_flags "${CMAKE_CXX_FLAGS}")
@@ -47,6 +58,8 @@ set(CMAKE_CXX_FLAGS "${_save_cxx_flags}")
 unset(_save_c_flags)
 unset(_save_cxx_flags)
 
+endif(USE_EXTERNAL_GTEST)
+
 #
 # Build PROJ unit tests
 #
diff --git a/test/unit/Makefile.am b/test/unit/Makefile.am
index 34624b0d..210bdc90 100644
--- a/test/unit/Makefile.am
+++ b/test/unit/Makefile.am
@@ -4,7 +4,7 @@ EXTRA_DIST = CMakeLists.txt
 
 noinst_HEADERS = gtest_include.h
 
-AM_CPPFLAGS = -I$(top_srcdir)/src -I$(top_srcdir)/include -I$(top_srcdir)/test -I$(top_srcdir)/test/googletest/include @SQLITE3_FLAGS@
+AM_CPPFLAGS = -I$(top_srcdir)/src -I$(top_srcdir)/include -I$(top_srcdir)/test @GTEST_CFLAGS@ @SQLITE3_CFLAGS@
 AM_CXXFLAGS = @CXX_WFLAGS@ @NO_ZERO_AS_NULL_POINTER_CONSTANT_FLAG@
 
 PROJ_LIB ?= ../../data
@@ -19,43 +19,43 @@ noinst_PROGRAMS += gie_self_tests
 noinst_PROGRAMS += include_proj_h_from_c
 
 pj_transform_test_SOURCES = pj_transform_test.cpp main.cpp
-pj_transform_test_LDADD = ../../src/libproj.la ../../test/googletest/libgtest.la
+pj_transform_test_LDADD = ../../src/libproj.la @GTEST_LIBS@
 
 pj_transform_test-check: pj_transform_test
 	PROJ_LIB=$(PROJ_LIB) ./pj_transform_test
 
 pj_phi2_test_SOURCES = pj_phi2_test.cpp main.cpp
-pj_phi2_test_LDADD = ../../src/libproj.la ../../test/googletest/libgtest.la
+pj_phi2_test_LDADD = ../../src/libproj.la @GTEST_LIBS@
 
 pj_phi2_test-check: pj_phi2_test
 	./pj_phi2_test
 
 proj_errno_string_test_SOURCES = proj_errno_string_test.cpp main.cpp
-proj_errno_string_test_LDADD= ../../src/libproj.la ../../test/googletest/libgtest.la
+proj_errno_string_test_LDADD= ../../src/libproj.la @GTEST_LIBS@
 
 proj_errno_string_test-check: proj_errno_string_test
 	./proj_errno_string_test
 
 proj_angular_io_test_SOURCES = proj_angular_io_test.cpp main.cpp
-proj_angular_io_test_LDADD = ../../src/libproj.la ../../test/googletest/libgtest.la
+proj_angular_io_test_LDADD = ../../src/libproj.la @GTEST_LIBS@
 
 proj_angular_io_test-check: proj_angular_io_test
 	./proj_angular_io_test
 
 proj_context_test_SOURCES = proj_context_test.cpp main.cpp
-proj_context_test_LDADD = ../../src/libproj.la ../../test/googletest/libgtest.la
+proj_context_test_LDADD = ../../src/libproj.la @GTEST_LIBS@
 
 proj_context_test-check: proj_context_test
 	./proj_context_test
 
 test_cpp_api_SOURCES = test_util.cpp test_common.cpp test_crs.cpp test_metadata.cpp test_io.cpp test_operation.cpp test_datum.cpp test_factory.cpp test_c_api.cpp main.cpp
-test_cpp_api_LDADD = ../../src/libproj.la ../../test/googletest/libgtest.la @SQLITE3_LDFLAGS@
+test_cpp_api_LDADD = ../../src/libproj.la @GTEST_LIBS@ @SQLITE3_LIBS@
 
 test_cpp_api-check: test_cpp_api
 	PROJ_LIB=$(PROJ_LIB) ./test_cpp_api
 
 gie_self_tests_SOURCES = gie_self_tests.cpp main.cpp
-gie_self_tests_LDADD = ../../src/libproj.la ../../test/googletest/libgtest.la @SQLITE3_LDFLAGS@
+gie_self_tests_LDADD = ../../src/libproj.la @GTEST_LIBS@ @SQLITE3_LIBS@
 
 gie_self_tests-check: gie_self_tests
 	PROJ_LIB=$(PROJ_LIB) ./gie_self_tests
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index ff29185f..d535e412 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -507,36 +507,21 @@ TEST_F(CApi, proj_as_proj_string_incompatible_WKT1) {
 
 // ---------------------------------------------------------------------------
 
-TEST_F(CApi, proj_as_proj_string_etmerc_option_yes) {
+TEST_F(CApi, proj_as_proj_string_approx_tmerc_option_yes) {
     auto obj = proj_create(m_ctxt, "+proj=tmerc +type=crs");
     ObjectKeeper keeper(obj);
     ASSERT_NE(obj, nullptr);
 
-    const char *options[] = {"USE_ETMERC=YES", nullptr};
+    const char *options[] = {"USE_APPROX_TMERC=YES", nullptr};
     auto str = proj_as_proj_string(m_ctxt, obj, PJ_PROJ_4, options);
     ASSERT_NE(str, nullptr);
     EXPECT_EQ(str,
-              std::string("+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 "
+              std::string("+proj=tmerc +approx +lat_0=0 +lon_0=0 +k=1 +x_0=0 "
                           "+y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"));
 }
 
 // ---------------------------------------------------------------------------
 
-TEST_F(CApi, proj_as_proj_string_etmerc_option_no) {
-    auto obj = proj_create(m_ctxt, "+proj=utm +zone=31 +type=crs");
-    ObjectKeeper keeper(obj);
-    ASSERT_NE(obj, nullptr);
-
-    const char *options[] = {"USE_ETMERC=NO", nullptr};
-    auto str = proj_as_proj_string(m_ctxt, obj, PJ_PROJ_4, options);
-    ASSERT_NE(str, nullptr);
-    EXPECT_EQ(str, std::string("+proj=tmerc +lat_0=0 +lon_0=3 +k=0.9996 "
-                               "+x_0=500000 +y_0=0 +datum=WGS84 +units=m "
-                               "+no_defs +type=crs"));
-}
-
-// ---------------------------------------------------------------------------
-
 TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) {
     auto crs = proj_create_from_database(m_ctxt, "EPSG", "3844",
                                          PJ_CATEGORY_CRS, false, nullptr);
@@ -3022,4 +3007,210 @@ TEST_F(CApi, proj_create_cartesian_2D_cs) {
     }
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_get_crs_info_list_from_database) {
+    { proj_crs_info_list_destroy(nullptr); }
+
+    { proj_get_crs_list_parameters_destroy(nullptr); }
+
+    // All null parameters
+    {
+        auto list = proj_get_crs_info_list_from_database(nullptr, nullptr,
+                                                         nullptr, nullptr);
+        ASSERT_NE(list, nullptr);
+        ASSERT_NE(list[0], nullptr);
+        EXPECT_NE(list[0]->auth_name, nullptr);
+        EXPECT_NE(list[0]->code, nullptr);
+        EXPECT_NE(list[0]->name, nullptr);
+        proj_crs_info_list_destroy(list);
+    }
+
+    // Default parameters
+    {
+        int result_count = 0;
+        auto params = proj_get_crs_list_parameters_create();
+        auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params,
+                                                         &result_count);
+        proj_get_crs_list_parameters_destroy(params);
+        ASSERT_NE(list, nullptr);
+        EXPECT_GT(result_count, 1);
+        EXPECT_EQ(list[result_count], nullptr);
+        bool found4326 = false;
+        bool found4978 = false;
+        bool found4979 = false;
+        bool found32631 = false;
+        bool found3855 = false;
+        bool found3901 = false;
+        for (int i = 0; i < result_count; i++) {
+            auto code = std::string(list[i]->code);
+            if (code == "4326") {
+                found4326 = true;
+                EXPECT_EQ(std::string(list[i]->auth_name), "EPSG");
+                EXPECT_EQ(std::string(list[i]->name), "WGS 84");
+                EXPECT_EQ(list[i]->type, PJ_TYPE_GEOGRAPHIC_2D_CRS);
+                EXPECT_EQ(list[i]->deprecated, 0);
+                EXPECT_EQ(list[i]->bbox_valid, 1);
+                EXPECT_EQ(list[i]->west_lon_degree, -180.0);
+                EXPECT_EQ(list[i]->south_lat_degree, -90.0);
+                EXPECT_EQ(list[i]->east_lon_degree, 180.0);
+                EXPECT_EQ(list[i]->north_lat_degree, 90.0);
+                EXPECT_EQ(std::string(list[i]->area_name), "World");
+                EXPECT_EQ(list[i]->projection_method_name, nullptr);
+            } else if (code == "4978") {
+                found4978 = true;
+                EXPECT_EQ(list[i]->type, PJ_TYPE_GEOCENTRIC_CRS);
+            } else if (code == "4979") {
+                found4979 = true;
+                EXPECT_EQ(list[i]->type, PJ_TYPE_GEOGRAPHIC_3D_CRS);
+            } else if (code == "32631") {
+                found32631 = true;
+                EXPECT_EQ(list[i]->type, PJ_TYPE_PROJECTED_CRS);
+                EXPECT_EQ(std::string(list[i]->projection_method_name),
+                          "Transverse Mercator");
+            } else if (code == "3855") {
+                found3855 = true;
+                EXPECT_EQ(list[i]->type, PJ_TYPE_VERTICAL_CRS);
+            } else if (code == "3901") {
+                found3901 = true;
+                EXPECT_EQ(list[i]->type, PJ_TYPE_COMPOUND_CRS);
+            }
+            EXPECT_EQ(list[i]->deprecated, 0);
+        }
+        EXPECT_TRUE(found4326);
+        EXPECT_TRUE(found4978);
+        EXPECT_TRUE(found4979);
+        EXPECT_TRUE(found32631);
+        EXPECT_TRUE(found3855);
+        EXPECT_TRUE(found3901);
+        proj_crs_info_list_destroy(list);
+    }
+
+    // Filter on only geodetic crs
+    {
+        int result_count = 0;
+        auto params = proj_get_crs_list_parameters_create();
+        params->typesCount = 1;
+        auto type = PJ_TYPE_GEODETIC_CRS;
+        params->types = &type;
+        auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params,
+                                                         &result_count);
+        bool foundGeog2D = false;
+        bool foundGeog3D = false;
+        bool foundGeocentric = false;
+        for (int i = 0; i < result_count; i++) {
+            foundGeog2D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS;
+            foundGeog3D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_3D_CRS;
+            foundGeocentric |= list[i]->type == PJ_TYPE_GEOCENTRIC_CRS;
+            EXPECT_TRUE(list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
+                        list[i]->type == PJ_TYPE_GEOGRAPHIC_3D_CRS ||
+                        list[i]->type == PJ_TYPE_GEOCENTRIC_CRS);
+        }
+        EXPECT_TRUE(foundGeog2D);
+        EXPECT_TRUE(foundGeog3D);
+        EXPECT_TRUE(foundGeocentric);
+        proj_get_crs_list_parameters_destroy(params);
+        proj_crs_info_list_destroy(list);
+    }
+
+    // Filter on only geographic crs
+    {
+        int result_count = 0;
+        auto params = proj_get_crs_list_parameters_create();
+        params->typesCount = 1;
+        auto type = PJ_TYPE_GEOGRAPHIC_CRS;
+        params->types = &type;
+        auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params,
+                                                         &result_count);
+        bool foundGeog2D = false;
+        bool foundGeog3D = false;
+        for (int i = 0; i < result_count; i++) {
+            foundGeog2D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS;
+            foundGeog3D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_3D_CRS;
+            EXPECT_TRUE(list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
+                        list[i]->type == PJ_TYPE_GEOGRAPHIC_3D_CRS);
+        }
+        EXPECT_TRUE(foundGeog2D);
+        EXPECT_TRUE(foundGeog3D);
+        proj_get_crs_list_parameters_destroy(params);
+        proj_crs_info_list_destroy(list);
+    }
+
+    // Filter on only geographic 2D crs and projected CRS
+    {
+        int result_count = 0;
+        auto params = proj_get_crs_list_parameters_create();
+        params->typesCount = 2;
+        const PJ_TYPE types[] = {PJ_TYPE_GEOGRAPHIC_2D_CRS,
+                                 PJ_TYPE_PROJECTED_CRS};
+        params->types = types;
+        auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params,
+                                                         &result_count);
+        bool foundGeog2D = false;
+        bool foundProjected = false;
+        for (int i = 0; i < result_count; i++) {
+            foundGeog2D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS;
+            foundProjected |= list[i]->type == PJ_TYPE_PROJECTED_CRS;
+            EXPECT_TRUE(list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
+                        list[i]->type == PJ_TYPE_PROJECTED_CRS);
+        }
+        EXPECT_TRUE(foundGeog2D);
+        EXPECT_TRUE(foundProjected);
+        proj_get_crs_list_parameters_destroy(params);
+        proj_crs_info_list_destroy(list);
+    }
+
+    // Filter on bbox (inclusion)
+    {
+        int result_count = 0;
+        auto params = proj_get_crs_list_parameters_create();
+        params->bbox_valid = 1;
+        params->west_lon_degree = 2;
+        params->south_lat_degree = 49;
+        params->east_lon_degree = 2.1;
+        params->north_lat_degree = 49.1;
+        params->typesCount = 1;
+        auto type = PJ_TYPE_PROJECTED_CRS;
+        params->types = &type;
+        auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params,
+                                                         &result_count);
+        ASSERT_NE(list, nullptr);
+        EXPECT_GT(result_count, 1);
+        for (int i = 0; i < result_count; i++) {
+            EXPECT_LE(list[i]->west_lon_degree, params->west_lon_degree);
+            EXPECT_LE(list[i]->south_lat_degree, params->south_lat_degree);
+            EXPECT_GE(list[i]->east_lon_degree, params->east_lon_degree);
+            EXPECT_GE(list[i]->north_lat_degree, params->north_lat_degree);
+        }
+        proj_get_crs_list_parameters_destroy(params);
+        proj_crs_info_list_destroy(list);
+    }
+
+    // Filter on bbox (intersection)
+    {
+        int result_count = 0;
+        auto params = proj_get_crs_list_parameters_create();
+        params->bbox_valid = 1;
+        params->west_lon_degree = 2;
+        params->south_lat_degree = 49;
+        params->east_lon_degree = 2.1;
+        params->north_lat_degree = 49.1;
+        params->crs_area_of_use_contains_bbox = 0;
+        params->typesCount = 1;
+        auto type = PJ_TYPE_PROJECTED_CRS;
+        params->types = &type;
+        auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params,
+                                                         &result_count);
+        ASSERT_NE(list, nullptr);
+        EXPECT_GT(result_count, 1);
+        for (int i = 0; i < result_count; i++) {
+            EXPECT_LE(list[i]->west_lon_degree, params->east_lon_degree);
+            EXPECT_LE(list[i]->south_lat_degree, params->north_lat_degree);
+            EXPECT_GE(list[i]->east_lon_degree, params->west_lon_degree);
+            EXPECT_GE(list[i]->north_lat_degree, params->south_lat_degree);
+        }
+        proj_get_crs_list_parameters_destroy(params);
+        proj_crs_info_list_destroy(list);
+    }
+}
 } // namespace
diff --git a/test/unit/test_datum.cpp b/test/unit/test_datum.cpp
index 18cf244a..c6ccf5cf 100644
--- a/test/unit/test_datum.cpp
+++ b/test/unit/test_datum.cpp
@@ -127,6 +127,18 @@ TEST(datum, ellipsoid_from_inverse_flattening) {
 
 // ---------------------------------------------------------------------------
 
+TEST(datum, ellipsoid_from_null_inverse_flattening) {
+
+    auto ellipsoid = Ellipsoid::createFlattenedSphere(
+        PropertyMap(), Length(6378137), Scale(0));
+    EXPECT_FALSE(ellipsoid->inverseFlattening().has_value());
+    EXPECT_FALSE(ellipsoid->semiMinorAxis().has_value());
+    EXPECT_FALSE(ellipsoid->semiMedianAxis().has_value());
+    EXPECT_TRUE(ellipsoid->isSphere());
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(datum, ellipsoid_from_semi_minor_axis) {
 
     auto ellipsoid = Ellipsoid::createTwoAxis(PropertyMap(), Length(6378137),
diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp
index 81f745f5..ce019079 100644
--- a/test/unit/test_factory.cpp
+++ b/test/unit/test_factory.cpp
@@ -222,6 +222,16 @@ TEST(factory, AuthorityFactory_createExtent) {
 
 // ---------------------------------------------------------------------------
 
+TEST(factory, AuthorityFactory_createExtent_no_bbox) {
+    auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto extent = factory->createExtent("1361"); // Sudan - south. Deprecated
+    EXPECT_EQ(*(extent->description()), "Sudan - south");
+    const auto &geogElts = extent->geographicElements();
+    EXPECT_TRUE(geogElts.empty());
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(factory, AuthorityFactory_createGeodeticDatum) {
     auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
     EXPECT_THROW(factory->createGeodeticDatum("-1"),
@@ -414,6 +424,17 @@ TEST(factory, AuthorityFactory_createGeodeticCRS_geographic2D) {
 
 // ---------------------------------------------------------------------------
 
+TEST(factory, AuthorityFactory_createGeodeticCRS_geographic2D_area_no_bbox) {
+    auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto crs = factory->createGeodeticCRS("4296"); // Sudan - deprecated
+    auto domain = crs->domains()[0];
+    auto extent = domain->domainOfValidity();
+    ASSERT_TRUE(extent != nullptr);
+    EXPECT_TRUE(extent->isEquivalentTo(factory->createExtent("1361").get()));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(factory, AuthorityFactory_createGeodeticCRS_geographic3D) {
     auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
     auto crs = factory->createGeodeticCRS("4979");
@@ -2789,4 +2810,103 @@ TEST(factory, listAreaOfUseFromName) {
     }
 }
 
+// ---------------------------------------------------------------------------
+
+TEST(factory, getCRSInfoList) {
+    auto ctxt = DatabaseContext::create();
+    {
+        auto factory = AuthorityFactory::create(ctxt, std::string());
+        auto list = factory->getCRSInfoList();
+        EXPECT_GT(list.size(), 1U);
+        bool foundEPSG = false;
+        bool foundIGNF = true;
+        bool found4326 = false;
+        for (const auto &info : list) {
+            foundEPSG |= info.authName == "EPSG";
+            foundIGNF |= info.authName == "IGNF";
+            if (info.authName == "EPSG" && info.code == "4326") {
+                found4326 = true;
+            }
+        }
+        EXPECT_TRUE(foundEPSG);
+        EXPECT_TRUE(foundIGNF);
+        EXPECT_TRUE(found4326);
+    }
+    {
+        auto factory = AuthorityFactory::create(ctxt, "EPSG");
+        auto list = factory->getCRSInfoList();
+        EXPECT_GT(list.size(), 1U);
+        bool found4326 = false;
+        bool found4978 = false;
+        bool found4979 = false;
+        bool found32631 = false;
+        bool found3855 = false;
+        bool found6871 = false;
+        for (const auto &info : list) {
+            EXPECT_EQ(info.authName, "EPSG");
+            if (info.code == "4326") {
+                EXPECT_EQ(info.name, "WGS 84");
+                EXPECT_EQ(info.type,
+                          AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS);
+                EXPECT_EQ(info.deprecated, false);
+                EXPECT_EQ(info.bbox_valid, true);
+                EXPECT_EQ(info.west_lon_degree, -180.0);
+                EXPECT_EQ(info.south_lat_degree, -90.0);
+                EXPECT_EQ(info.east_lon_degree, 180.0);
+                EXPECT_EQ(info.north_lat_degree, 90.0);
+                EXPECT_EQ(info.areaName, "World");
+                EXPECT_TRUE(info.projectionMethodName.empty());
+                found4326 = true;
+            } else if (info.code == "4296") { // Soudan - deprecated
+                EXPECT_EQ(info.bbox_valid, false);
+                EXPECT_EQ(info.west_lon_degree, 0.0);
+                EXPECT_EQ(info.south_lat_degree, 0.0);
+                EXPECT_EQ(info.east_lon_degree, 0.0);
+                EXPECT_EQ(info.north_lat_degree, 0.0);
+            } else if (info.code == "4978") {
+                EXPECT_EQ(info.name, "WGS 84");
+                EXPECT_EQ(info.type,
+                          AuthorityFactory::ObjectType::GEOCENTRIC_CRS);
+                found4978 = true;
+            } else if (info.code == "4979") {
+                EXPECT_EQ(info.name, "WGS 84");
+                EXPECT_EQ(info.type,
+                          AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS);
+                found4979 = true;
+            } else if (info.code == "32631") {
+                EXPECT_EQ(info.name, "WGS 84 / UTM zone 31N");
+                EXPECT_EQ(info.type,
+                          AuthorityFactory::ObjectType::PROJECTED_CRS);
+                EXPECT_EQ(info.deprecated, false);
+                EXPECT_EQ(info.bbox_valid, true);
+                EXPECT_EQ(info.west_lon_degree, 0.0);
+                EXPECT_EQ(info.south_lat_degree, 0.0);
+                EXPECT_EQ(info.east_lon_degree, 6.0);
+                EXPECT_EQ(info.north_lat_degree, 84.0);
+                EXPECT_EQ(info.areaName, "World - N hemisphere - 0\xC2\xB0"
+                                         "E to 6\xC2\xB0"
+                                         "E - by country");
+                EXPECT_EQ(info.projectionMethodName, "Transverse Mercator");
+                found32631 = true;
+            } else if (info.code == "3855") {
+                EXPECT_EQ(info.name, "EGM2008 height");
+                EXPECT_EQ(info.type,
+                          AuthorityFactory::ObjectType::VERTICAL_CRS);
+                found3855 = true;
+            } else if (info.code == "6871") {
+                EXPECT_EQ(info.name,
+                          "WGS 84 / Pseudo-Mercator +  EGM2008 geoid height");
+                EXPECT_EQ(info.type,
+                          AuthorityFactory::ObjectType::COMPOUND_CRS);
+                found6871 = true;
+            }
+        }
+        EXPECT_TRUE(found4326);
+        EXPECT_TRUE(found4978);
+        EXPECT_TRUE(found4979);
+        EXPECT_TRUE(found32631);
+        EXPECT_TRUE(found3855);
+        EXPECT_TRUE(found6871);
+    }
+}
 } // namespace
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 1bbedc03..3cb50352 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -7654,12 +7654,39 @@ TEST(io, projparse_etmerc) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 "
+        "+proj=tmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 "
         "+datum=WGS84 +units=m +no_defs +type=crs");
 
     auto wkt1 = crs->exportToWKT(
         WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
-    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=etmerc +lat_0=0 "
+    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\"") == std::string::npos)
+        << wkt1;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_tmerc_approx) {
+    auto obj =
+        PROJStringParser().createFromPROJString("+proj=tmerc +approx +type=crs");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto wkt2 = crs->exportToWKT(
+        &WKTFormatter::create()->simulCurNodeHasId().setMultiLine(false));
+    EXPECT_TRUE(
+        wkt2.find("METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]]") !=
+        std::string::npos)
+        << wkt2;
+
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=tmerc +approx +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 "
+        "+datum=WGS84 +units=m +no_defs +type=crs");
+
+    auto wkt1 = crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=tmerc +approx +lat_0=0 "
                           "+lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m "
                           "+no_defs\"]") != std::string::npos)
         << wkt1;
@@ -8735,6 +8762,10 @@ TEST(io, projparse_errors) {
     EXPECT_THROW(PROJStringParser().createFromPROJString(
                      "proj=pipeline step init=epsg:4326 init=epsg:4326"),
                  ParsingException);
+
+    EXPECT_THROW(
+        PROJStringParser().createFromPROJString("proj=\tinit= +type=crs"),
+        ParsingException);
 }
 
 // ---------------------------------------------------------------------------
@@ -8961,3 +8992,61 @@ TEST(wkt_export, precision) {
             WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
         wkt);
 }
+
+// ---------------------------------------------------------------------------
+
+// Avoid division by zero
+
+TEST(wkt_export, invalid_linear_unit) {
+    auto wkt = "PROJCS[\"WGS 84 / UTM zone 31N\",\n"
+               "    GEOGCS[\"WGS 84\",\n"
+               "        DATUM[\"WGS_1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563]],\n"
+               "        PRIMEM[\"Greenwich\",0],\n"
+               "        UNIT[\"degree\",0.0174532925199433]],\n"
+               "    PROJECTION[\"Transverse_Mercator\"],\n"
+               "    PARAMETER[\"latitude_of_origin\",0],\n"
+               "    PARAMETER[\"central_meridian\",3],\n"
+               "    PARAMETER[\"scale_factor\",0.9996],\n"
+               "    PARAMETER[\"false_easting\",500000],\n"
+               "    PARAMETER[\"false_northing\",0],\n"
+               "    UNIT[\"foo\",0]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    EXPECT_THROW(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
+        FormattingException);
+}
+
+// ---------------------------------------------------------------------------
+
+// Avoid division by zero
+
+TEST(wkt_export, invalid_angular_unit) {
+    auto wkt = "PROJCS[\"WGS 84 / UTM zone 31N\",\n"
+               "    GEOGCS[\"WGS 84\",\n"
+               "        DATUM[\"WGS_1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563]],\n"
+               "        PRIMEM[\"Greenwich\",0],\n"
+               "        UNIT[\"foo\",0]],\n"
+               "    PROJECTION[\"Transverse_Mercator\"],\n"
+               "    PARAMETER[\"latitude_of_origin\",0],\n"
+               "    PARAMETER[\"central_meridian\",3],\n"
+               "    PARAMETER[\"scale_factor\",0.9996],\n"
+               "    PARAMETER[\"false_easting\",500000],\n"
+               "    PARAMETER[\"false_northing\",0],\n"
+               "    UNIT[\"meter\",1]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    EXPECT_THROW(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
+        FormattingException);
+}
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 6c1ecf1c..586226e2 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -1237,9 +1237,9 @@ TEST(operation, tmerc_export) {
 
     {
         auto formatter = PROJStringFormatter::create();
-        formatter->setUseETMercForTMerc(true);
+        formatter->setUseApproxTMerc(true);
         EXPECT_EQ(conv->exportToPROJString(formatter.get()),
-                  "+proj=etmerc +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5");
+                  "+proj=tmerc +approx +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5");
     }
 
     EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()),
@@ -6061,21 +6061,21 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) {
               "+ellps=WGS84");
     {
         auto formatter = PROJStringFormatter::create();
-        formatter->setUseETMercForTMerc(true);
+        formatter->setUseApproxTMerc(true);
         EXPECT_EQ(op->exportToPROJString(formatter.get()),
-                  "+proj=pipeline +step +inv +proj=etmerc +lat_0=1 +lon_0=2 "
+                  "+proj=pipeline +step +inv +proj=tmerc +approx +lat_0=1 +lon_0=2 "
                   "+k=3 +x_0=4 +y_0=5 +ellps=WGS84 +step "
                   "+proj=vgridshift +grids=bla.gtx +multiplier=0.001 +step "
-                  "+proj=utm +zone=32 "
+                  "+proj=utm +approx +zone=32 "
                   "+ellps=WGS84");
     }
     {
         auto formatter = PROJStringFormatter::create();
-        formatter->setUseETMercForTMerc(true);
+        formatter->setUseApproxTMerc(true);
         EXPECT_EQ(op->inverse()->exportToPROJString(formatter.get()),
-                  "+proj=pipeline +step +inv +proj=utm +zone=32 +ellps=WGS84 "
+                  "+proj=pipeline +step +inv +proj=utm +approx +zone=32 +ellps=WGS84 "
                   "+step +inv +proj=vgridshift +grids=bla.gtx "
-                  "+multiplier=0.001 +step +proj=etmerc +lat_0=1 +lon_0=2 "
+                  "+multiplier=0.001 +step +proj=tmerc +approx +lat_0=1 +lon_0=2 "
                   "+k=3 +x_0=4 +y_0=5 +ellps=WGS84");
     }
 
diff --git a/travis/mingw32/install.sh b/travis/mingw32/install.sh
index defe1485..e4fcfe72 100755
--- a/travis/mingw32/install.sh
+++ b/travis/mingw32/install.sh
@@ -13,15 +13,15 @@ ln -s /usr/x86_64-w64-mingw32/lib/libwinpthread-1.dll  $HOME/.wine/drive_c/windo
 # build sqlite3
 wget https://sqlite.org/2018/sqlite-autoconf-3250100.tar.gz
 tar xzf sqlite-autoconf-3250100.tar.gz
-(cd sqlite-autoconf-3250100 &&  CC="ccache x86_64-w64-mingw32-gcc" CXX="ccache x86_64-w64-mingw32-g++" LD=x86_64-w64-mingw32-ld ./configure --host=x86_64-w64-mingw32 --prefix=/tmp/proj_autoconf_install && make -j3 && make install)
-ln -s /tmp/proj_autoconf_install/bin/libsqlite3-0.dll  $HOME/.wine/drive_c/windows
+(cd sqlite-autoconf-3250100 &&  CC="ccache x86_64-w64-mingw32-gcc" CXX="ccache x86_64-w64-mingw32-g++" LD=x86_64-w64-mingw32-ld ./configure --host=x86_64-w64-mingw32 --prefix=/usr/x86_64-w64-mingw32 && make -j3 && sudo make install)
+ln -s /usr/x86_64-w64-mingw32/bin/libsqlite3-0.dll  $HOME/.wine/drive_c/windows
 
 # prepare build files
 ./autogen.sh
 # autoconf build
 mkdir build_autoconf
 cd build_autoconf
-CC="ccache x86_64-w64-mingw32-gcc" CXX="ccache x86_64-w64-mingw32-g++" LD=x86_64-w64-mingw32-ld ../configure --host=x86_64-w64-mingw32 --with-sqlite3_include=/tmp/proj_autoconf_install/include --with-sqlite3_ldflags="-L/tmp/proj_autoconf_install/lib -lsqlite3" --prefix=/tmp/proj_autoconf_install
+CC="ccache x86_64-w64-mingw32-gcc" CXX="ccache x86_64-w64-mingw32-g++" LD=x86_64-w64-mingw32-ld ../configure --host=x86_64-w64-mingw32 --prefix=/tmp/proj_autoconf_install
 make -j3
 make install
 make dist-all
@@ -37,7 +37,7 @@ cd ..
 # autoconf build with grids
 mkdir build_autoconf_grids
 cd build_autoconf_grids
-CC="ccache x86_64-w64-mingw32-gcc" CXX="ccache x86_64-w64-mingw32-g++" LD=x86_64-w64-mingw32-ld ../configure --host=x86_64-w64-mingw32 --with-sqlite3_include=/tmp/proj_autoconf_install/include --with-sqlite3_ldflags="-L/tmp/proj_autoconf_install/lib -lsqlite3" --prefix=/tmp/proj_autoconf_install_nad
+CC="ccache x86_64-w64-mingw32-gcc" CXX="ccache x86_64-w64-mingw32-g++" LD=x86_64-w64-mingw32-ld ../configure --host=x86_64-w64-mingw32 --prefix=/tmp/proj_autoconf_install_nad
 make -j3
 make install
 find /tmp/proj_autoconf_install_nad