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

27 files changed, 977 insertions, 1461 deletions

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