Merge pull request #1855 from rouault/rfc4_merge_back_master

[RFC4_dev] Merge master back to rfc4 latest branch

22 files changed, 713 insertions, 110 deletions

diff --git a/test/cli/CMakeLists.txt b/test/cli/CMakeLists.txt
index 4e1ab75a..feef5bf0 100644
--- a/test/cli/CMakeLists.txt
+++ b/test/cli/CMakeLists.txt
@@ -7,6 +7,7 @@ set(PROJINFO_BIN "projinfo")
 set(CCT_BIN "cct")
 proj_add_test_script_sh("test27" PROJ_BIN)
 proj_add_test_script_sh("test83" PROJ_BIN)
+proj_add_test_script_sh("testproj" PROJ_BIN)
 proj_add_test_script_sh("testvarious" CS2CS_BIN)
 proj_add_test_script_sh("testdatumfile" CS2CS_BIN "connu")
 proj_add_test_script_sh("testIGNF" CS2CS_BIN "ntf_r93.gsb")
diff --git a/test/cli/Makefile.am b/test/cli/Makefile.am
index 253b85d8..42ee466d 100644
--- a/test/cli/Makefile.am
+++ b/test/cli/Makefile.am
@@ -1,5 +1,5 @@
 # Executables paths passed to test scripts
-PROJ_LIB ?= ../../data
+PROJ_LIB = ../../data/for_tests
 THIS_DIR = $(top_srcdir)/test/cli
 EXEPATH = ../../src
 PROJEXE = $(EXEPATH)/proj
@@ -10,6 +10,7 @@ CCTEXE = $(EXEPATH)/cct
 # PROJ.4 test scripts
 TEST27 = $(THIS_DIR)/test27
 TEST83 = $(THIS_DIR)/test83
+TESTPROJ = $(THIS_DIR)/testproj
 TESTNTV2 = $(THIS_DIR)/testntv2
 TESTVARIOUS = $(THIS_DIR)/testvarious
 TESTFLAKY = $(THIS_DIR)/testflaky
@@ -24,6 +25,7 @@ EXTRA_DIST = pj_out27.dist pj_out83.dist td_out.dist \
 		testIGNF proj_outIGNF.dist \
 		testprojinfo testprojinfo_out.dist \
 		testcct testcct_out.dist \
+		testproj testproj_out.dist \
 		CMakeLists.txt
 
 testprojinfo-check:
@@ -35,6 +37,9 @@ test27-check:
 test83-check:
 	$(TEST83) $(PROJEXE)
 
+testproj-check:
+	$(TESTPROJ) $(PROJEXE)
+
 testvarious-check:
 	PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(TESTVARIOUS) $(CS2CSEXE)
 
@@ -56,4 +61,4 @@ testntv2-check:
 testcct-check:
 	PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(TESTCCT) $(CCTEXE)
 
-check-local: testprojinfo-check test27-check test83-check testvarious-check testdatumfile-check testign-check testntv2-check testcct-check
+check-local: testprojinfo-check test27-check test83-check testproj-check testvarious-check testdatumfile-check testign-check testntv2-check testcct-check
diff --git a/test/cli/ntv2_out.dist b/test/cli/ntv2_out.dist
index ce866f37..d90fdc37 100644
--- a/test/cli/ntv2_out.dist
+++ b/test/cli/ntv2_out.dist
@@ -20,3 +20,10 @@ Interpolating very close (and sometimes a bit outside) to the edges a NTv2 subgr
 -115.49166667 51.07500000 0	-115.49062909	51.07497666 0.00000000
 -115.49166668 51.07500000 0	-115.49062910	51.07497666 0.00000000
 -115.49166667 51.07499990 0	-115.49062909	51.07497656 0.00000000
+##############################################################
+Attempt first with ntv2_0.gsb and then conus
+-111.5 45.26	-111.50079772	45.25992835 0.00000000
+##############################################################
+NAD27 -> NAD83: 1st through ntv2, 2nd through conus
+55d00'00.000"N 111d00'00.000"W 0.0	55d0'0.366"N	111d0'3.251"W 0.000
+39d00'00.000"N 111d00'00.000"W 0.0	38d59'59.912"N	111d0'2.604"W 0.000
diff --git a/test/cli/td_out.dist b/test/cli/td_out.dist
index 76c6e6ec..cf4b8d73 100644
--- a/test/cli/td_out.dist
+++ b/test/cli/td_out.dist
@@ -24,9 +24,9 @@ edge or even a wee bit outside (#141).
 -5.4999 51.9999	-5.498793593803	52.000009531513 0.000000000000
 -5.5001 52.0	-5.500100000000	52.000000000000 0.000000000000
 ##############################################################
-NAD27 -> NAD83: 1st through ntv1, 2nd through conus
-55d00'00.000"N 111d00'00.000"W 0.0	55d0'0.339"N	111d0'3.245"W 0.000
-39d00'00.000"N 111d00'00.000"W 0.0	38d59'59.912"N	111d0'2.604"W 0.000
+NAD27 -> NAD83: 1st through ntv1 or ntv2, 2nd through conus
+55d00'00.000"N 111d00'00.000"W 0.0	55.0001	-111.0009 0.0000
+39d00'00.000"N 111d00'00.000"W 0.0	39.0000	-111.0007 0.0000
 ##############################################################
 WGS84 -> WGS84+EGM96
 2dE 49dN 0	2.00	49.00 -44.64
diff --git a/test/cli/testdatumfile b/test/cli/testdatumfile
index 9bd12ce4..16e4bbc3 100755
--- a/test/cli/testdatumfile
+++ b/test/cli/testdatumfile
@@ -102,9 +102,9 @@ $EXE +proj=latlong +datum=WGS84 \
 EOF
 #
 echo "##############################################################" >> ${OUT}
-echo "NAD27 -> NAD83: 1st through ntv1, 2nd through conus" >> ${OUT}
+echo "NAD27 -> NAD83: 1st through ntv1 or ntv2, 2nd through conus" >> ${OUT}
 #
-$EXE NAD27 NAD83 -E >>${OUT} <<EOF
+$EXE NAD27 NAD83 -E -f "%.4f" >>${OUT} <<EOF
 55d00'00.000"N 111d00'00.000"W 0.0
 39d00'00.000"N 111d00'00.000"W 0.0
 EOF
diff --git a/test/cli/testntv2 b/test/cli/testntv2
index d15ad7aa..44ccac1e 100755
--- a/test/cli/testntv2
+++ b/test/cli/testntv2
@@ -71,6 +71,20 @@ $EXE +proj=latlong +datum=NAD83 +to +proj=latlong +ellps=clrk66 +nadgrids=ntv2_0
 -115.49166667 51.07499990 0
 EOF
 
+echo "##############################################################" >> ${OUT}
+echo Attempt first with ntv2_0.gsb and then conus >> ${OUT}
+$EXE +proj=longlat +datum=NAD27 +to +proj=longlat +datum=WGS84 -E -d 8 >>${OUT} <<EOF
+-111.5 45.26
+EOF
+
+echo "##############################################################" >> ${OUT}
+echo "NAD27 -> NAD83: 1st through ntv2, 2nd through conus" >> ${OUT}
+#
+$EXE NAD27 NAD83 -E >>${OUT} <<EOF
+55d00'00.000"N 111d00'00.000"W 0.0
+39d00'00.000"N 111d00'00.000"W 0.0
+EOF
+
 #
 ##############################################################################
 # Done! 
diff --git a/test/cli/testproj b/test/cli/testproj
new file mode 100755
index 00000000..8686224e
--- /dev/null
+++ b/test/cli/testproj
@@ -0,0 +1,55 @@
+:
+# Script to test proj exe
+#
+TEST_CLI_DIR=`dirname $0`
+EXE=$1
+
+usage()
+{
+    echo "Usage: ${0} <path to 'proj' program>" 
+    echo
+    exit 1
+}
+
+if test -z "${EXE}"; then
+    EXE=../../src/cs2cs
+fi
+
+if test ! -x ${EXE}; then
+    echo "*** ERROR: Can not find '${EXE}' program!"
+    exit 1
+fi
+
+if test -z "${PROJ_LIB}"; then
+    export PROJ_LIB="`dirname $0`/../../data"
+fi
+
+echo "============================================"
+echo "Running ${0} using ${EXE}:"
+echo "============================================"
+
+OUT=testproj_out
+#
+echo "doing tests into file ${OUT}, please wait"
+#
+$EXE +ellps=WGS84 +proj=ob_tran +o_proj=latlon +o_lon_p=0.0 +o_lat_p=90.0 +lon_0=360.0 +to_meter=0.0174532925199433 +no_defs -E -f '%.3f' >${OUT} <<EOF
+2 49
+EOF
+
+#
+# do 'diff' with distribution results
+echo "diff ${OUT} with testproj_out.dist"
+diff -u -b ${OUT} ${TEST_CLI_DIR}/testproj_out.dist
+if [ $? -ne 0 ] ; then
+	echo  ""
+	echo "PROBLEMS HAVE OCCURRED"
+	echo "test file ${OUT} saved"
+    echo
+	exit 100
+else
+	echo "TEST OK"
+	echo "test file ${OUT} removed"
+    echo
+	/bin/rm -f ${OUT}
+	exit 0
+fi
diff --git a/test/cli/testproj_out.dist b/test/cli/testproj_out.dist
new file mode 100644
index 00000000..fad0a014
--- /dev/null
+++ b/test/cli/testproj_out.dist
@@ -0,0 +1 @@
+2 49	2.000	49.000
diff --git a/test/cli/testprojinfo b/test/cli/testprojinfo
index 4ce5e90c..b9c452fb 100755
--- a/test/cli/testprojinfo
+++ b/test/cli/testprojinfo
@@ -155,6 +155,26 @@ echo 'Testing -s "GDA2020" -t "AHD height" --grid-check none -o PROJ --spatial-t
 $EXE -s "GDA2020" -t "AHD height" --grid-check none -o PROJ --spatial-test intersects >>${OUT} 2>&1
 echo "" >>${OUT}
 
+echo 'Testing -k ellipsoid WGS84' >> ${OUT}
+$EXE -k ellipsoid WGS84 >>${OUT} 2>&1
+echo "" >>${OUT}
+
+echo 'Testing -k ellipsoid EPSG:7030' >> ${OUT}
+$EXE -k ellipsoid EPSG:7030 >>${OUT} 2>&1
+echo "" >>${OUT}
+
+echo 'Testing -k datum WGS84' >> ${OUT}
+$EXE -k datum WGS84 >>${OUT} 2>&1
+echo "" >>${OUT}
+
+echo 'Testing -k datum EPSG:6326' >> ${OUT}
+$EXE -k datum EPSG:6326 >>${OUT} 2>&1
+echo "" >>${OUT}
+
+echo 'Testing -k operation EPSG:8457 -o PROJ -q' >> ${OUT}
+$EXE -k operation EPSG:8457 -o PROJ -q >>${OUT} 2>&1
+echo "" >>${OUT}
+
 # do 'diff' with distribution results
 echo "diff ${OUT} with testprojinfo_out.dist"
 diff -u ${OUT} ${TEST_CLI_DIR}/testprojinfo_out.dist
diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
index 25b29a79..03ce9046 100644
--- a/test/cli/testprojinfo_out.dist
+++ b/test/cli/testprojinfo_out.dist
@@ -182,7 +182,7 @@ PROJJSON:
 Testing projinfo -s EPSG:4326 -t EPSG:32631
 Candidate operations found: 1
 -------------------------------------
-Operation n°1:
+Operation No. 1:
 
 EPSG:16031, UTM zone 31N, 0 m, World - N hemisphere - 0°E to 6°E
 
@@ -214,7 +214,7 @@ Testing projinfo -s NAD27 -t NAD83
 Candidate operations found: 1
 Note: using '--spatial-test intersects' would bring more results (10)
 -------------------------------------
-Operation n°1:
+Operation No. 1:
 
 unknown id, Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, World, has ballpark transformation
 
@@ -268,12 +268,12 @@ COORDINATEOPERATION["Ballpark geographic offset from NAD27 to NAD83",
 
 Testing projinfo -s NAD27 -t NAD83 --grid-check none --spatial-test intersects --summary
 Candidate operations found: 10
-DERIVED_FROM(EPSG):1312, NAD27 to NAD83 (3), 1.0 m, Canada
+DERIVED_FROM(EPSG):1312, NAD27 to NAD83 (3), 2.0 m, Canada
 DERIVED_FROM(EPSG):1313, NAD27 to NAD83 (4), 1.5 m, Canada - NAD27
 DERIVED_FROM(EPSG):1241, NAD27 to NAD83 (1), 0.15 m, USA - CONUS including EEZ
 DERIVED_FROM(EPSG):1243, NAD27 to NAD83 (2), 0.5 m, USA - Alaska including EEZ
-EPSG:1462, NAD27 to NAD83 (5), 1.0 m, Canada - Quebec
-EPSG:1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
+DERIVED_FROM(EPSG):1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
+EPSG:1462, NAD27 to NAD83 (5), 2.0 m, Canada - Quebec
 EPSG:9111, NAD27 to NAD83 (9), 1.5 m, Canada - Saskatchewan
 unknown id, Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, World, has ballpark transformation
 EPSG:8555, NAD27 to NAD83 (7), 0.15 m, USA - CONUS and GoM
@@ -282,9 +282,9 @@ EPSG:8549, NAD27 to NAD83 (8), 0.5 m, USA - Alaska
 Testing projinfo -s NAD27 -t NAD83 --grid-check none --spatial-test intersects
 Candidate operations found: 10
 -------------------------------------
-Operation n°1:
+Operation No. 1:
 
-DERIVED_FROM(EPSG):1312, NAD27 to NAD83 (3), 1.0 m, Canada
+DERIVED_FROM(EPSG):1312, NAD27 to NAD83 (3), 2.0 m, Canada
 
 PROJ string:
 +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=ntv1_can.dat +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
@@ -324,7 +324,7 @@ COORDINATEOPERATION["NAD27 to NAD83 (3)",
     METHOD["NTv1",
         ID["EPSG",9614]],
     PARAMETERFILE["Latitude and longitude difference file","ntv1_can.dat"],
-    OPERATIONACCURACY[1.0],
+    OPERATIONACCURACY[2.0],
     USAGE[
         SCOPE["Historic record only - now superseded - see remarks."],
         AREA["Canada"],
@@ -332,7 +332,7 @@ COORDINATEOPERATION["NAD27 to NAD83 (3)",
     ID["DERIVED_FROM(EPSG)",1312]]
 
 -------------------------------------
-Operation n°2:
+Operation No. 2:
 
 DERIVED_FROM(EPSG):1313, NAD27 to NAD83 (4), 1.5 m, Canada - NAD27
 
@@ -382,7 +382,7 @@ COORDINATEOPERATION["NAD27 to NAD83 (4)",
     ID["DERIVED_FROM(EPSG)",1313]]
 
 -------------------------------------
-Operation n°3:
+Operation No. 3:
 
 DERIVED_FROM(EPSG):1241, NAD27 to NAD83 (1), 0.15 m, USA - CONUS including EEZ
 
@@ -431,7 +431,7 @@ COORDINATEOPERATION["NAD27 to NAD83 (1)",
     ID["DERIVED_FROM(EPSG)",1241]]
 
 -------------------------------------
-Operation n°4:
+Operation No. 4:
 
 DERIVED_FROM(EPSG):1243, NAD27 to NAD83 (2), 0.5 m, USA - Alaska including EEZ
 
@@ -480,16 +480,15 @@ COORDINATEOPERATION["NAD27 to NAD83 (2)",
     ID["DERIVED_FROM(EPSG)",1243]]
 
 -------------------------------------
-Operation n°5:
+Operation No. 5:
 
-EPSG:1462, NAD27 to NAD83 (5), 1.0 m, Canada - Quebec
+DERIVED_FROM(EPSG):1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
 
 PROJ string:
-+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GS2783v1.QUE +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=na27na83.gsb +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 WKT2:2019 string:
-COORDINATEOPERATION["NAD27 to NAD83 (5)",
-    VERSION["SGQ-Can QC NT1"],
+COORDINATEOPERATION["NAD27 to NAD83 (6)",
     SOURCECRS[
         GEOGCRS["NAD27",
             DATUM["North American Datum 1927",
@@ -520,28 +519,27 @@ COORDINATEOPERATION["NAD27 to NAD83 (5)",
                     ORDER[2],
                     ANGLEUNIT["degree",0.0174532925199433]],
             ID["EPSG",4269]]],
-    METHOD["NTv1",
-        ID["EPSG",9614]],
-    PARAMETERFILE["Latitude and longitude difference file","GS2783v1.QUE"],
-    OPERATIONACCURACY[1.0],
+    METHOD["NTv2",
+        ID["EPSG",9615]],
+    PARAMETERFILE["Latitude and longitude difference file","na27na83.gsb"],
+    OPERATIONACCURACY[1.5],
     USAGE[
-        SCOPE["Historic record only - now superseded - see remarks."],
+        SCOPE["Accuracy 1-2 metres."],
         AREA["Canada - Quebec"],
         BBOX[44.99,-79.85,62.62,-57.1]],
-    ID["EPSG",1462],
-    REMARK["Densification for Quebec of code 1312. Replaced by NAD27 to NAD83 (6) (code 1573). Uses NT method which expects longitudes positive west; EPSG GeogCRSs NAD27 (code 4267) and NAD83 (code 4269) have longitudes positive east."]]
+    ID["DERIVED_FROM(EPSG)",1573]]
 
 -------------------------------------
-Operation n°6:
+Operation No. 6:
 
-EPSG:1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
+EPSG:1462, NAD27 to NAD83 (5), 2.0 m, Canada - Quebec
 
 PROJ string:
-+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=NA27NA83.GSB +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GS2783v1.QUE +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 WKT2:2019 string:
-COORDINATEOPERATION["NAD27 to NAD83 (6)",
-    VERSION["SGQ-Can QC NT2"],
+COORDINATEOPERATION["NAD27 to NAD83 (5)",
+    VERSION["SGQ-Can QC NT1"],
     SOURCECRS[
         GEOGCRS["NAD27",
             DATUM["North American Datum 1927",
@@ -572,19 +570,19 @@ COORDINATEOPERATION["NAD27 to NAD83 (6)",
                     ORDER[2],
                     ANGLEUNIT["degree",0.0174532925199433]],
             ID["EPSG",4269]]],
-    METHOD["NTv2",
-        ID["EPSG",9615]],
-    PARAMETERFILE["Latitude and longitude difference file","NA27NA83.GSB"],
-    OPERATIONACCURACY[1.5],
+    METHOD["NTv1",
+        ID["EPSG",9614]],
+    PARAMETERFILE["Latitude and longitude difference file","GS2783v1.QUE"],
+    OPERATIONACCURACY[2.0],
     USAGE[
-        SCOPE["Accuracy 1-2 metres."],
+        SCOPE["Historic record only - now superseded - see remarks."],
         AREA["Canada - Quebec"],
         BBOX[44.99,-79.85,62.62,-57.1]],
-    ID["EPSG",1573],
-    REMARK["Also distributed with file name QUE27-83.gsb. Replaces NAD27 to NAD83 (5) (code 1462). Uses NT method which expects longitudes positive west; EPSG GeogCRSs NAD27 (code 4267) and NAD83 (code 4269) have longitudes positive east."]]
+    ID["EPSG",1462],
+    REMARK["Densification for Quebec of code 1312. Replaced by NAD27 to NAD83 (6) (code 1573). Uses NT method which expects longitudes positive west; EPSG GeogCRSs NAD27 (code 4267) and NAD83 (code 4269) have longitudes positive east."]]
 
 -------------------------------------
-Operation n°7:
+Operation No. 7:
 
 EPSG:9111, NAD27 to NAD83 (9), 1.5 m, Canada - Saskatchewan
 
@@ -635,7 +633,7 @@ COORDINATEOPERATION["NAD27 to NAD83 (9)",
     ID["EPSG",9111]]
 
 -------------------------------------
-Operation n°8:
+Operation No. 8:
 
 unknown id, Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, World, has ballpark transformation
 
@@ -688,7 +686,7 @@ COORDINATEOPERATION["Ballpark geographic offset from NAD27 to NAD83",
         BBOX[-90,-180,90,180]]]
 
 -------------------------------------
-Operation n°9:
+Operation No. 9:
 
 EPSG:8555, NAD27 to NAD83 (7), 0.15 m, USA - CONUS and GoM
 
@@ -740,7 +738,7 @@ COORDINATEOPERATION["NAD27 to NAD83 (7)",
     REMARK["Uses NADCON5 method which expects longitudes positive east in range 0-360°; source and target CRSs have longitudes positive east in range -180° to +180°. Accuracy at 67% confidence level is 0.15m onshore, 1m nearshore and undetermined farther offshore."]]
 
 -------------------------------------
-Operation n°10:
+Operation No. 10:
 
 EPSG:8549, NAD27 to NAD83 (8), 0.5 m, USA - Alaska
 
@@ -903,7 +901,7 @@ PROJ.4 string:
 Testing RH2000 height to SWEREF99: projinfo -s EPSG:5613 -t EPSG:4977
 Candidate operations found: 1
 -------------------------------------
-Operation n°1:
+Operation No. 1:
 
 PROJ:EPSG_5613_TO_EPSG_4977, RH2000 height to SWEREF99, unknown accuracy, Sweden - onshore
 
@@ -949,7 +947,7 @@ COORDINATEOPERATION["RH2000 height to SWEREF99",
 Testing NAD83(2011) + NAVD88 height -> NAD83(2011) : projinfo -s EPSG:6349 -t EPSG:6319 --spatial-test intersects -o PROJ
 Candidate operations found: 2
 -------------------------------------
-Operation n°1:
+Operation No. 1:
 
 unknown id, Inverse of NAD83(2011) to NAVD88 height (3), 0.015 m, USA - CONUS - onshore
 
@@ -957,7 +955,7 @@ PROJ string:
 +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=g2018u0.gtx +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 -------------------------------------
-Operation n°2:
+Operation No. 2:
 
 unknown id, Inverse of NAD83(2011) to NAVD88 height (2), 0.02 m, USA - Alaska
 
@@ -967,7 +965,7 @@ PROJ string:
 Testing NGF IGN69 height to RGF93: projinfo -s EPSG:5720 -t EPSG:4965 -o PROJ
 Candidate operations found: 2
 -------------------------------------
-Operation n°1:
+Operation No. 1:
 
 INVERSE(DERIVED_FROM(EPSG)):8885, Inverse of RGF93 to NGF IGN69 height (3), 0.01 m, France - mainland onshore
 
@@ -975,7 +973,7 @@ PROJ string:
 +proj=vgridshift +grids=RAF18.gtx +multiplier=1
 
 -------------------------------------
-Operation n°2:
+Operation No. 2:
 
 INVERSE(EPSG):10000, Inverse of RGF93 to NGF IGN69 height (1), 0.5 m, France - mainland onshore
 
@@ -1081,7 +1079,7 @@ unknown id, Inverse of 3-degree Gauss-Kruger zone 3 + DHDN to ETRS89 (2), 3 m, G
 Testing -s "GDA94" -t "AHD height" --grid-check none -o PROJ --spatial-test intersects
 Candidate operations found: 1
 -------------------------------------
-Operation n°1:
+Operation No. 1:
 
 DERIVED_FROM(EPSG):5656, GDA94 to AHD height (49), 0.03 m, Australia - mainland
 
@@ -1091,10 +1089,45 @@ PROJ string:
 Testing -s "GDA2020" -t "AHD height" --grid-check none -o PROJ --spatial-test intersects
 Candidate operations found: 1
 -------------------------------------
-Operation n°1:
+Operation No. 1:
 
 DERIVED_FROM(EPSG):8451, GDA2020 to AHD height (1), 0.03 m, Australia Christmas and Cocos - onshore
 
 PROJ string:
 +proj=pipeline +step +inv +proj=vgridshift +grids=AUSGeoid2020_20180201.gtx +multiplier=1
 
+Testing -k ellipsoid WGS84
+PROJ string:
++ellps=WGS84
+
+WKT2:2019 string:
+ELLIPSOID["WGS 84",6378137,298.257223563,
+    LENGTHUNIT["metre",1],
+    ID["EPSG",7030]]
+
+Testing -k ellipsoid EPSG:7030
+PROJ string:
++ellps=WGS84
+
+WKT2:2019 string:
+ELLIPSOID["WGS 84",6378137,298.257223563,
+    LENGTHUNIT["metre",1],
+    ID["EPSG",7030]]
+
+Testing -k datum WGS84
+WKT2:2019 string:
+DATUM["World Geodetic System 1984",
+    ELLIPSOID["WGS 84",6378137,298.257223563,
+        LENGTHUNIT["metre",1]],
+    ID["EPSG",6326]]
+
+Testing -k datum EPSG:6326
+WKT2:2019 string:
+DATUM["World Geodetic System 1984",
+    ELLIPSOID["WGS 84",6378137,298.257223563,
+        LENGTHUNIT["metre",1]],
+    ID["EPSG",6326]]
+
+Testing -k operation EPSG:8457 -o PROJ -q
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=bessel +step +proj=helmert +x=674.374 +y=15.056 +z=405.346 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
+
diff --git a/test/fuzzers/standard_fuzzer.cpp b/test/fuzzers/standard_fuzzer.cpp
index 29ea50f8..468e8cbb 100644
--- a/test/fuzzers/standard_fuzzer.cpp
+++ b/test/fuzzers/standard_fuzzer.cpp
@@ -34,7 +34,8 @@
 #include <sys/types.h>
 #include <unistd.h>
 
-#include "proj_internal.h" // For pj_gc_unloadall()
+#define ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
+#include "proj.h"
 #include "proj_api.h"
 
 /* Standalone build:
@@ -105,10 +106,7 @@ int LLVMFuzzerTestOneInput(const uint8_t *buf, size_t len)
     {
         free(buf_dup);
         pj_free(pj_src);
-#ifndef OMIT_DEALLOCATION
-        pj_gc_unloadall(pj_get_default_ctx());
-        pj_deallocate_grids();
-#endif
+        proj_cleanup();
         return 0;
     }
     double x = 0, y = 0, z = 9;
@@ -135,10 +133,7 @@ int LLVMFuzzerTestOneInput(const uint8_t *buf, size_t len)
         free(buf_dup);
         pj_free(pj_src);
         pj_free(pj_dst);
-#ifndef OMIT_DEALLOCATION
-        pj_gc_unloadall(pj_get_default_ctx());
-        pj_deallocate_grids();
-#endif
+        proj_cleanup();
         return 0;
     }
 #ifdef STANDALONE
@@ -161,10 +156,7 @@ int LLVMFuzzerTestOneInput(const uint8_t *buf, size_t len)
     free(buf_dup);
     pj_free(pj_src);
     pj_free(pj_dst);
-#ifndef OMIT_DEALLOCATION
-    pj_gc_unloadall(pj_get_default_ctx());
-    pj_deallocate_grids();
-#endif
+    proj_cleanup();
     return 0;
 }
 
diff --git a/test/gie/Makefile.am b/test/gie/Makefile.am
index 1539cb2e..783dfb7c 100644
--- a/test/gie/Makefile.am
+++ b/test/gie/Makefile.am
@@ -12,7 +12,7 @@ EXTRA_DIST = 4D-API_cs2cs-style.gie \
 	DHDN_ETRS89.gie \
 	geotiff_grids.gie
 
-PROJ_LIB ?= ../../data
+PROJ_LIB ?= ../../data/for_tests
 
 4D-API-cs2cs-style: 4D-API_cs2cs-style.gie
 	PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(GIEEXE) $<
diff --git a/test/gigs/Makefile.am b/test/gigs/Makefile.am
index 56493550..22ab099b 100644
--- a/test/gigs/Makefile.am
+++ b/test/gigs/Makefile.am
@@ -20,7 +20,7 @@ EXTRA_DIST = \
 	5201.gie \
 	5208.gie
 
-PROJ_LIB ?= ../../data
+PROJ_LIB ?= ../../data/for_tests
 
 5101.1: 5101.1-jhs.gie
 	PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(GIEEXE) $<
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index e14d4c70..f5890f1a 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -71,7 +71,7 @@ target_link_libraries(proj_pj_transform_test
   ${PROJ_LIBRARIES})
 add_test(NAME proj_pj_transform_test COMMAND proj_pj_transform_test)
 set_property(TEST proj_pj_transform_test
-        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data")
+        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data/for_tests")
 
 
 add_executable(proj_errno_string_test
@@ -82,7 +82,7 @@ target_link_libraries(proj_errno_string_test
   ${PROJ_LIBRARIES})
 add_test(NAME proj_errno_string_test COMMAND proj_errno_string_test)
 set_property(TEST proj_errno_string_test
-        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data")
+        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data/for_tests")
 
 add_executable(proj_angular_io_test
   main.cpp
@@ -92,7 +92,7 @@ target_link_libraries(proj_angular_io_test
   ${PROJ_LIBRARIES})
 add_test(NAME proj_angular_io_test COMMAND proj_angular_io_test)
 set_property(TEST proj_angular_io_test
-        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data")
+        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data/for_tests")
 
 add_executable(proj_context_test
   main.cpp
@@ -102,7 +102,7 @@ target_link_libraries(proj_context_test
   ${PROJ_LIBRARIES})
 add_test(NAME proj_context_test COMMAND proj_context_test)
 set_property(TEST proj_context_test
-        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data")
+        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data/for_tests")
 
 if(MSVC AND BUILD_LIBPROJ_SHARED)
   # ph_phi2_test not compatible of a .dll build
@@ -115,7 +115,7 @@ else()
     ${PROJ_LIBRARIES})
   add_test(NAME pj_phi2_test COMMAND pj_phi2_test)
   set_property(TEST pj_phi2_test
-        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data")
+        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data/for_tests")
 endif()
 
 add_executable(proj_test_cpp_api
@@ -135,13 +135,8 @@ target_link_libraries(proj_test_cpp_api
   ${PROJ_LIBRARIES}
   ${SQLITE3_LIBRARY})
 add_test(NAME proj_test_cpp_api COMMAND proj_test_cpp_api)
-if(MSVC)
-    set_property(TEST proj_test_cpp_api
-            PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data\\;${PROJECT_SOURCE_DIR}/data")
-else()
-    set_property(TEST proj_test_cpp_api
-            PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data:${PROJECT_SOURCE_DIR}/data")
-endif()
+set_property(TEST proj_test_cpp_api
+        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data/for_tests")
 
 
 add_executable(gie_self_tests
@@ -152,7 +147,7 @@ target_link_libraries(gie_self_tests
   ${PROJ_LIBRARIES})
 add_test(NAME gie_self_tests COMMAND gie_self_tests)
 set_property(TEST gie_self_tests
-        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data")
+        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data/for_tests")
 
 
 add_executable(test_network
@@ -167,10 +162,5 @@ target_link_libraries(test_network
   ${PROJ_LIBRARIES}
   ${SQLITE3_LIBRARY})
 add_test(NAME test_network COMMAND test_network)
-if(MSVC)
-    set_property(TEST test_network
-        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data\\;${PROJECT_SOURCE_DIR}/data;PROJ_SOURCE_DATA=${PROJECT_SOURCE_DIR}/data")
-else()
-    set_property(TEST test_network
-        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data:${PROJECT_SOURCE_DIR}/data;PROJ_SOURCE_DATA=${PROJECT_SOURCE_DIR}/data")
-endif()
+set_property(TEST test_network
+        PROPERTY ENVIRONMENT "PROJ_IGNORE_USER_WRITABLE_DIRECTORY=YES;PROJ_LIB=${PROJECT_BINARY_DIR}/data/for_tests;PROJ_SOURCE_DATA=${PROJECT_SOURCE_DIR}/data")
diff --git a/test/unit/Makefile.am b/test/unit/Makefile.am
index b7a26d4e..86ed383a 100644
--- a/test/unit/Makefile.am
+++ b/test/unit/Makefile.am
@@ -7,7 +7,7 @@ noinst_HEADERS = gtest_include.h test_primitives.hpp
 AM_CPPFLAGS = -I$(top_srcdir)/src -I$(top_srcdir)/include -I$(top_srcdir)/test @GTEST_CFLAGS@ @SQLITE3_CFLAGS@
 AM_CXXFLAGS = @CXX_WFLAGS@ @NO_ZERO_AS_NULL_POINTER_CONSTANT_FLAG@
 
-PROJ_LIB ?= ../../data
+PROJ_LIB ?= ../../data/for_tests
 
 noinst_PROGRAMS = pj_transform_test
 noinst_PROGRAMS += pj_phi2_test
diff --git a/test/unit/pj_transform_test.cpp b/test/unit/pj_transform_test.cpp
index 5ca8dcce..ddb054f0 100644
--- a/test/unit/pj_transform_test.cpp
+++ b/test/unit/pj_transform_test.cpp
@@ -632,6 +632,23 @@ TEST(proj_api_h, default_fileapi) {
 
 // ---------------------------------------------------------------------------
 
+TEST(pj_transform_test, ob_tran_to_meter_as_dest) {
+    auto src = pj_init_plus(
+        "+ellps=WGS84 +a=57.29577951308232 +proj=eqc +lon_0=0.0 +no_defs");
+    auto dst = pj_init_plus("+ellps=WGS84 +proj=ob_tran +o_proj=latlon "
+                            "+o_lon_p=0.0 +o_lat_p=90.0 +lon_0=360.0 "
+                            "+to_meter=0.0174532925199433 +no_defs");
+    double x = 2 * DEG_TO_RAD;
+    double y = 49 * DEG_TO_RAD;
+    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
+    EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
+    EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
+    pj_free(src);
+    pj_free(dst);
+}
+
+// ---------------------------------------------------------------------------
+
 struct Spy {
     bool gotInMyFOpen = false;
     bool gotInMyFRead = false;
@@ -705,4 +722,19 @@ TEST(proj_api_h, custom_fileapi) {
     EXPECT_TRUE(spy.gotInMyFClose);
 }
 
+TEST(pj_transform_test, ob_tran_to_meter_as_srouce) {
+    auto src = pj_init_plus("+ellps=WGS84 +proj=ob_tran +o_proj=latlon "
+                            "+o_lon_p=0.0 +o_lat_p=90.0 +lon_0=360.0 "
+                            "+to_meter=0.0174532925199433 +no_defs");
+    auto dst = pj_init_plus(
+        "+ellps=WGS84 +a=57.29577951308232 +proj=eqc +lon_0=0.0 +no_defs");
+    double x = 2 * DEG_TO_RAD;
+    double y = 49 * DEG_TO_RAD;
+    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
+    EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
+    EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
+    pj_free(src);
+    pj_free(dst);
+}
+
 } // namespace
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 8871f679..5c6dabba 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -1668,7 +1668,7 @@ TEST_F(CApi, proj_create_from_name) {
                                          false, 0, nullptr);
         ASSERT_NE(res, nullptr);
         ObjListKeeper keeper_res(res);
-        EXPECT_EQ(proj_list_get_count(res), 4);
+        EXPECT_EQ(proj_list_get_count(res), 5);
     }
     {
         auto res = proj_create_from_name(m_ctxt, "xx", "WGS 84", nullptr, 0,
@@ -3486,6 +3486,36 @@ TEST_F(CApi, proj_normalize_for_visualization_on_crs) {
 
 // ---------------------------------------------------------------------------
 
+TEST_F(CApi, proj_coordoperation_create_inverse) {
+
+    auto P = proj_create(
+        m_ctxt, "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+                "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push "
+                "+v_3 +step +proj=cart +ellps=evrst30 +step +proj=helmert "
+                "+x=293 +y=836 +z=318 +rx=0.5 +ry=1.6 +rz=-2.8 +s=2.1 "
+                "+convention=position_vector +step +inv +proj=cart "
+                "+ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert "
+                "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto Pinversed = proj_coordoperation_create_inverse(m_ctxt, P);
+    ObjectKeeper keeper_Pinversed(Pinversed);
+    ASSERT_NE(Pinversed, nullptr);
+
+    auto projstr = proj_as_proj_string(m_ctxt, Pinversed, PJ_PROJ_5, nullptr);
+    ASSERT_NE(projstr, nullptr);
+    EXPECT_EQ(std::string(projstr),
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=293 "
+              "+y=836 +z=318 +rx=0.5 +ry=1.6 +rz=-2.8 +s=2.1 "
+              "+convention=position_vector +step +inv +proj=cart "
+              "+ellps=evrst30 +step +proj=pop +v_3 +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST_F(CApi, proj_get_remarks) {
     auto co = proj_create_from_database(m_ctxt, "EPSG", "8048",
                                         PJ_CATEGORY_COORDINATE_OPERATION, false,
@@ -4456,4 +4486,34 @@ TEST_F(CApi, proj_context_set_sqlite3_vfs_name) {
     proj_context_destroy(ctx);
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_is_equivalent_to_with_ctx) {
+    auto from_epsg = proj_create_from_database(m_ctxt, "EPSG", "7844",
+                                               PJ_CATEGORY_CRS, false, nullptr);
+    ObjectKeeper keeper_from_epsg(from_epsg);
+    ASSERT_NE(from_epsg, nullptr);
+
+    auto wkt = "GEOGCRS[\"GDA2020\",\n"
+               "    DATUM[\"GDA2020\",\n"
+               "        ELLIPSOID[\"GRS_1980\",6378137,298.257222101,\n"
+               "            LENGTHUNIT[\"metre\",1]]],\n"
+               "    PRIMEM[\"Greenwich\",0,\n"
+               "        ANGLEUNIT[\"Degree\",0.0174532925199433]],\n"
+               "    CS[ellipsoidal,2],\n"
+               "        AXIS[\"geodetic latitude (Lat)\",north,\n"
+               "            ORDER[1],\n"
+               "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+               "        AXIS[\"geodetic longitude (Lon)\",east,\n"
+               "            ORDER[2],\n"
+               "            ANGLEUNIT[\"degree\",0.0174532925199433]]]";
+    auto from_wkt =
+        proj_create_from_wkt(m_ctxt, wkt, nullptr, nullptr, nullptr);
+    ObjectKeeper keeper_from_wkt(from_wkt);
+    EXPECT_NE(from_wkt, nullptr);
+
+    EXPECT_TRUE(proj_is_equivalent_to_with_ctx(m_ctxt, from_epsg, from_wkt,
+                                               PJ_COMP_EQUIVALENT));
+}
+
 } // namespace
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 427fd741..e6189a97 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -57,7 +57,8 @@ namespace {
 struct UnrelatedObject : public IComparable {
     UnrelatedObject() = default;
 
-    bool _isEquivalentTo(const IComparable *, Criterion) const override {
+    bool _isEquivalentTo(const IComparable *, Criterion,
+                         const DatabaseContextPtr &) const override {
         assert(false);
         return false;
     }
@@ -755,7 +756,9 @@ TEST(crs, EPSG_27561_projected_with_geodetic_in_grad_as_PROJ_string_and_WKT1) {
         nn_crs->DerivedCRS::isEquivalentTo(createUnrelatedObject().get()));
 
     auto wkt1 = crs->exportToWKT(
-        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL,
+                             DatabaseContext::create())
+            .get());
     EXPECT_EQ(
         wkt1,
         "PROJCS[\"NTF (Paris) / Lambert Nord France\",\n"
@@ -856,7 +859,9 @@ TEST(crs, EPSG_2222_projected_unit_foot_as_PROJ_string_and_WKT1) {
               "+x_0=213360 +y_0=0 +datum=NAD83 +units=ft +no_defs +type=crs");
 
     auto wkt1 = crs->exportToWKT(
-        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL,
+                             DatabaseContext::create())
+            .get());
     EXPECT_EQ(wkt1,
               "PROJCS[\"NAD83 / Arizona East (ft)\",\n"
               "    GEOGCS[\"NAD83\",\n"
@@ -903,11 +908,13 @@ TEST(crs, projected_with_parameter_unit_different_than_cs_unit_as_WKT1) {
     ASSERT_TRUE(crs != nullptr);
 
     auto wkt1 = crs->exportToWKT(
-        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL,
+                             DatabaseContext::create())
+            .get());
     EXPECT_EQ(wkt1,
               "PROJCS[\"unknown\",\n"
               "    GEOGCS[\"unknown\",\n"
-              "        DATUM[\"Unknown_based_on_GRS80_ellipsoid\",\n"
+              "        DATUM[\"Unknown based on GRS80 ellipsoid\",\n"
               "            SPHEROID[\"GRS 1980\",6378137,298.257222101]],\n"
               "        PRIMEM[\"Greenwich\",0],\n"
               "        UNIT[\"degree\",0.0174532925199433,\n"
@@ -943,6 +950,16 @@ TEST(crs, EPSG_32661_projected_north_pole_north_east) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               proj_string);
+
+    auto opNormalized = op->normalizeForVisualization();
+    auto proj_string_normalized =
+        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step "
+        "+proj=stere +lat_0=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 "
+        "+ellps=WGS84";
+    EXPECT_EQ(
+        opNormalized->exportToPROJString(PROJStringFormatter::create().get()),
+        proj_string_normalized);
 }
 
 // ---------------------------------------------------------------------------
@@ -964,6 +981,16 @@ TEST(crs, EPSG_5041_projected_north_pole_east_north) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               proj_string);
+
+    auto opNormalized = op->normalizeForVisualization();
+    auto proj_string_normalized =
+        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step "
+        "+proj=stere +lat_0=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 "
+        "+ellps=WGS84";
+    EXPECT_EQ(
+        opNormalized->exportToPROJString(PROJStringFormatter::create().get()),
+        proj_string_normalized);
 }
 
 // ---------------------------------------------------------------------------
@@ -985,6 +1012,16 @@ TEST(crs, EPSG_32761_projected_south_pole_north_east) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               proj_string);
+
+    auto opNormalized = op->normalizeForVisualization();
+    auto proj_string_normalized =
+        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step "
+        "+proj=stere +lat_0=-90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 "
+        "+ellps=WGS84";
+    EXPECT_EQ(
+        opNormalized->exportToPROJString(PROJStringFormatter::create().get()),
+        proj_string_normalized);
 }
 
 // ---------------------------------------------------------------------------
@@ -1006,6 +1043,16 @@ TEST(crs, EPSG_5042_projected_south_pole_east_north) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               proj_string);
+
+    auto opNormalized = op->normalizeForVisualization();
+    auto proj_string_normalized =
+        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step "
+        "+proj=stere +lat_0=-90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 "
+        "+ellps=WGS84";
+    EXPECT_EQ(
+        opNormalized->exportToPROJString(PROJStringFormatter::create().get()),
+        proj_string_normalized);
 }
 
 // ---------------------------------------------------------------------------
@@ -1511,6 +1558,75 @@ TEST(crs, geodeticcrs_identify_db) {
         EXPECT_EQ(res.front().first->identifiers()[0]->code(), "4326");
         EXPECT_EQ(res.front().second, 70);
     }
+
+    {
+        // Identify by code, but datum name is an alias of the official one
+        auto wkt = "GEOGCRS[\"GDA2020\",\n"
+                   "    DATUM[\"GDA2020\",\n"
+                   "        ELLIPSOID[\"GRS_1980\",6378137,298.257222101,\n"
+                   "            LENGTHUNIT[\"metre\",1]]],\n"
+                   "    PRIMEM[\"Greenwich\",0,\n"
+                   "        ANGLEUNIT[\"Degree\",0.0174532925199433]],\n"
+                   "    CS[ellipsoidal,2],\n"
+                   "        AXIS[\"geodetic latitude (Lat)\",north,\n"
+                   "            ORDER[1],\n"
+                   "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+                   "        AXIS[\"geodetic longitude (Lon)\",east,\n"
+                   "            ORDER[2],\n"
+                   "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+                   "    ID[\"EPSG\",7844]]";
+        auto obj = WKTParser().createFromWKT(wkt);
+        auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+
+        auto allFactory = AuthorityFactory::create(dbContext, std::string());
+        auto res = crs->identify(allFactory);
+        ASSERT_EQ(res.size(), 1U);
+        ASSERT_TRUE(!res.front().first->identifiers().empty());
+        EXPECT_EQ(*res.front().first->identifiers()[0]->codeSpace(), "EPSG");
+        EXPECT_EQ(res.front().first->identifiers()[0]->code(), "7844");
+        EXPECT_EQ(res.front().second, 100);
+
+        EXPECT_TRUE(crs->_isEquivalentTo(res.front().first.get(),
+                                         IComparable::Criterion::EQUIVALENT,
+                                         dbContext));
+        EXPECT_TRUE(res.front().first->_isEquivalentTo(
+            crs.get(), IComparable::Criterion::EQUIVALENT, dbContext));
+    }
+
+    {
+        // Identify by name, but datum name is an alias of the official one
+        auto wkt = "GEOGCRS[\"GDA2020\",\n"
+                   "    DATUM[\"GDA2020\",\n"
+                   "        ELLIPSOID[\"GRS_1980\",6378137,298.257222101,\n"
+                   "            LENGTHUNIT[\"metre\",1]]],\n"
+                   "    PRIMEM[\"Greenwich\",0,\n"
+                   "        ANGLEUNIT[\"Degree\",0.0174532925199433]],\n"
+                   "    CS[ellipsoidal,2],\n"
+                   "        AXIS[\"geodetic latitude (Lat)\",north,\n"
+                   "            ORDER[1],\n"
+                   "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+                   "        AXIS[\"geodetic longitude (Lon)\",east,\n"
+                   "            ORDER[2],\n"
+                   "            ANGLEUNIT[\"degree\",0.0174532925199433]]]";
+        auto obj = WKTParser().createFromWKT(wkt);
+        auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+
+        auto allFactory = AuthorityFactory::create(dbContext, std::string());
+        auto res = crs->identify(allFactory);
+        ASSERT_EQ(res.size(), 1U);
+        ASSERT_TRUE(!res.front().first->identifiers().empty());
+        EXPECT_EQ(*res.front().first->identifiers()[0]->codeSpace(), "EPSG");
+        EXPECT_EQ(res.front().first->identifiers()[0]->code(), "7844");
+        EXPECT_EQ(res.front().second, 100);
+
+        EXPECT_TRUE(crs->_isEquivalentTo(res.front().first.get(),
+                                         IComparable::Criterion::EQUIVALENT,
+                                         dbContext));
+        EXPECT_TRUE(res.front().first->_isEquivalentTo(
+            crs.get(), IComparable::Criterion::EQUIVALENT, dbContext));
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -2281,6 +2397,38 @@ TEST(crs, projectedCRS_identify_db) {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, projectedCRS_identify_wrong_auth_name_case) {
+    auto dbContext = DatabaseContext::create();
+    auto factoryAnonymous = AuthorityFactory::create(dbContext, std::string());
+    auto obj =
+        WKTParser()
+            .attachDatabaseContext(dbContext)
+            .setStrict(false)
+            .createFromWKT(
+                "PROJCS[\"World_Cylindrical_Equal_Area\","
+                "GEOGCS[\"GCS_WGS_1984\","
+                "DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\","
+                "6378137.0,298.257223563]],"
+                "PRIMEM[\"Greenwich\",0.0],"
+                "UNIT[\"Degree\",0.0174532925199433]],"
+                "PROJECTION[\"Cylindrical_Equal_Area\"],"
+                "PARAMETER[\"False_Easting\",0.0],"
+                "PARAMETER[\"False_Northing\",0.0],"
+                "PARAMETER[\"Central_Meridian\",0.0],"
+                "PARAMETER[\"Standard_Parallel_1\",0.0],UNIT[\"Meter\",1.0],"
+                "AUTHORITY[\"Esri\",54034]]"); // should be ESRI all caps
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto res = crs->identify(factoryAnonymous);
+    ASSERT_EQ(res.size(), 1U);
+    const auto &ids = res.front().first->identifiers();
+    ASSERT_EQ(ids.size(), 1U);
+    EXPECT_EQ(*(ids.front()->codeSpace()), "ESRI");
+    EXPECT_EQ(ids.front()->code(), "54034");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, mercator_1SP_as_WKT1_ESRI) {
 
     auto obj = PROJStringParser().createFromPROJString(
@@ -3351,6 +3499,17 @@ TEST(crs, compoundCRS_identify_db) {
         EXPECT_EQ(res.front().first->getEPSGCode(), 8769);
         EXPECT_EQ(res.front().second, 70);
     }
+    {
+        auto obj = PROJStringParser().createFromPROJString(
+            "+proj=tmerc +lat_0=0 +lon_0=72.05 +k=1 +x_0=3500000 "
+            "+y_0=-5811057.63 +ellps=krass "
+            "+towgs84=23.57,-140.95,-79.8,0,-0.35,-0.79,-0.22 "
+            "+geoidgrids=egm08_25.gtx +units=m +no_defs +type=crs");
+        auto crs = nn_dynamic_pointer_cast<CompoundCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        // Just check we don't get an exception
+        crs->identify(factory);
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -3753,6 +3912,21 @@ TEST(crs, boundCRS_identify_db) {
             WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
         EXPECT_TRUE(wkt.find("32122") != std::string::npos) << wkt;
     }
+
+    {
+        // Identify from a PROJ string with +towgs84
+        auto obj = PROJStringParser().createFromPROJString(
+            "+proj=utm +zone=48 +a=6377276.345 +b=6356075.41314024 "
+            "+towgs84=198,881,317,0,0,0,0 +units=m +no_defs +type=crs");
+        auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1U);
+        auto boundCRS = dynamic_cast<const BoundCRS *>(res.front().first.get());
+        ASSERT_TRUE(boundCRS != nullptr);
+        EXPECT_EQ(boundCRS->baseCRS()->getEPSGCode(), 3148);
+        EXPECT_EQ(res.front().second, 70);
+    }
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_datum.cpp b/test/unit/test_datum.cpp
index 8a5d67cf..adf0ae4b 100644
--- a/test/unit/test_datum.cpp
+++ b/test/unit/test_datum.cpp
@@ -44,7 +44,8 @@ namespace {
 struct UnrelatedObject : public IComparable {
     UnrelatedObject() = default;
 
-    bool _isEquivalentTo(const IComparable *, Criterion) const override {
+    bool _isEquivalentTo(const IComparable *, Criterion,
+                         const DatabaseContextPtr &) const override {
         assert(false);
         return false;
     }
diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp
index e342dad9..6a88e2b6 100644
--- a/test/unit/test_factory.cpp
+++ b/test/unit/test_factory.cpp
@@ -2737,8 +2737,8 @@ TEST(factory, createObjectsFromName) {
 
     EXPECT_EQ(factory->createObjectsFromName("").size(), 0U);
 
-    // ellipsoid + 3 geodeticCRS
-    EXPECT_EQ(factory->createObjectsFromName("WGS 84", {}, false).size(), 4U);
+    // ellipsoid + datum + 3 geodeticCRS
+    EXPECT_EQ(factory->createObjectsFromName("WGS 84", {}, false).size(), 5U);
 
     EXPECT_EQ(factory->createObjectsFromName("WGS 84", {}, true, 10).size(),
               10U);
@@ -2753,7 +2753,9 @@ TEST(factory, createObjectsFromName) {
         auto res = factoryEPSG->createObjectsFromName(
             "WGS84", {AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS}, true);
         EXPECT_EQ(res.size(),
-                  8U); // EPSG:4326 and EPSG:4030 and the 6 WGS84 realizations
+                  9U); // EPSG:4326 and EPSG:4030 and the 6 WGS84 realizations
+                       // and EPSG:7881 'Tritan St. Helena'' whose alias is
+                       // 'WGS 84 Tritan St. Helena'
         if (!res.empty()) {
             EXPECT_EQ(res.front()->getEPSGCode(), 4326);
         }
@@ -2876,7 +2878,7 @@ TEST(factory, getMetadata) {
     EXPECT_EQ(ctxt->getMetadata("i_do_not_exist"), nullptr);
     const char *IGNF_VERSION = ctxt->getMetadata("IGNF.VERSION");
     ASSERT_TRUE(IGNF_VERSION != nullptr);
-    EXPECT_EQ(std::string(IGNF_VERSION), "3.0.3");
+    EXPECT_EQ(std::string(IGNF_VERSION), "3.1.0");
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 38dfc2b4..fd38847c 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -493,18 +493,14 @@ TEST(wkt_parse, wkt1_geographic_old_datum_name_from_EPSG_code) {
 
 // ---------------------------------------------------------------------------
 
-TEST(wkt_parse, wkt1_geographic_old_datum_name_witout_EPSG_code) {
+TEST(wkt_parse, wkt1_geographic_old_datum_name_without_EPSG_code) {
     auto wkt =
         "GEOGCS[\"S-JTSK (Ferro)\",\n"
         "    "
         "DATUM[\"System_Jednotne_Trigonometricke_Site_Katastralni_Ferro\",\n"
-        "        SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,\n"
-        "            AUTHORITY[\"EPSG\",\"7004\"]]],\n"
-        "    PRIMEM[\"Ferro\",-17.66666666666667,\n"
-        "       AUTHORITY[\"EPSG\",\"8909\"]],\n"
-        "    UNIT[\"degree\",0.0174532925199433,\n"
-        "        AUTHORITY[\"EPSG\",\"9122\"]],\n"
-        "    AUTHORITY[\"EPSG\",\"4818\"]]";
+        "        SPHEROID[\"Bessel 1841\",6377397.155,299.1528128]],\n"
+        "    PRIMEM[\"Ferro\",-17.66666666666667],\n"
+        "    UNIT[\"degree\",0.0174532925199433]]";
     auto obj = WKTParser()
                    .attachDatabaseContext(DatabaseContext::create())
                    .createFromWKT(wkt);
@@ -5577,6 +5573,17 @@ TEST(wkt_parse, invalid_GEOCCS) {
                                            "NORTH],AXIS[\"longitude\",EAST]]"),
                  ParsingException);
 
+    // ellipsoidal CS is invalid in a GEOCCS
+    EXPECT_THROW(WKTParser().createFromWKT(
+                     "GEOCCS[\"WGS 84\",DATUM[\"World Geodetic System 1984\","
+                     "ELLIPSOID[\"WGS 84\",6378274,298.257223564,"
+                     "LENGTHUNIT[\"metre\",1]]],"
+                     "CS[ellipsoidal,2],AXIS[\"geodetic latitude (Lat)\",north,"
+                     "ANGLEUNIT[\"degree\",0.0174532925199433]],"
+                     "AXIS[\"geodetic longitude (Lon)\",east,"
+                     "ANGLEUNIT[\"degree\",0.0174532925199433]]]"),
+                 ParsingException);
+
     // 3 axis required
     EXPECT_THROW(WKTParser().createFromWKT(
                      "GEOCCS[\"x\",DATUM[\"x\",SPHEROID[\"x\",1,0.5]],PRIMEM["
@@ -6856,6 +6863,30 @@ TEST(io, projstringformatter_optim_hgridshift_vgridshift_hgridshift_inv) {
                   "+step +proj=pop +v_1 +v_2");
     }
 
+    // Test omit_fwd->omit_inv when inversing the pipeline
+    {
+        auto fmt = PROJStringFormatter::create();
+        fmt->startInversion();
+        fmt->ingestPROJString("+proj=hgridshift +grids=foo +omit_fwd");
+        fmt->stopInversion();
+
+        EXPECT_EQ(fmt->toString(),
+                  "+proj=pipeline "
+                  "+step +inv +proj=hgridshift +grids=foo +omit_inv");
+    }
+
+    // Test omit_inv->omit_fwd when inversing the pipeline
+    {
+        auto fmt = PROJStringFormatter::create();
+        fmt->startInversion();
+        fmt->ingestPROJString("+proj=hgridshift +grids=foo +omit_inv");
+        fmt->stopInversion();
+
+        EXPECT_EQ(fmt->toString(),
+                  "+proj=pipeline "
+                  "+step +inv +proj=hgridshift +grids=foo +omit_fwd");
+    }
+
     // Variant with first hgridshift inverted, and second forward
     {
         auto fmt = PROJStringFormatter::create();
@@ -9324,6 +9355,20 @@ TEST(io, createFromUserInput) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, createFromUserInput_hack_EPSG_102100) {
+    auto dbContext = DatabaseContext::create();
+    auto obj = createFromUserInput("EPSG:102100", dbContext);
+    auto crs = nn_dynamic_pointer_cast<CRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    const auto &ids = crs->identifiers();
+    ASSERT_EQ(ids.size(), 1U);
+    // we do not lie on the real authority
+    EXPECT_EQ(*ids[0]->codeSpace(), "ESRI");
+    EXPECT_EQ(ids[0]->code(), "102100");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, guessDialect) {
     EXPECT_EQ(WKTParser().guessDialect("LOCAL_CS[\"foo\"]"),
               WKTParser::WKTGuessedDialect::WKT1_GDAL);
@@ -11866,3 +11911,40 @@ TEST(json_import, multiple_ids) {
     EXPECT_EQ(ellps->exportToJSON(&(JSONFormatter::create()->setSchema("foo"))),
               json);
 }
+
+// ---------------------------------------------------------------------------
+
+TEST(json_export, coordinate_system_id) {
+    auto json = "{\n"
+                "  \"$schema\": \"foo\",\n"
+                "  \"type\": \"CoordinateSystem\",\n"
+                "  \"subtype\": \"ellipsoidal\",\n"
+                "  \"axis\": [\n"
+                "    {\n"
+                "      \"name\": \"Geodetic latitude\",\n"
+                "      \"abbreviation\": \"Lat\",\n"
+                "      \"direction\": \"north\",\n"
+                "      \"unit\": \"degree\"\n"
+                "    },\n"
+                "    {\n"
+                "      \"name\": \"Geodetic longitude\",\n"
+                "      \"abbreviation\": \"Lon\",\n"
+                "      \"direction\": \"east\",\n"
+                "      \"unit\": \"degree\"\n"
+                "    }\n"
+                "  ],\n"
+                "  \"id\": {\n"
+                "    \"authority\": \"EPSG\",\n"
+                "    \"code\": 6422\n"
+                "  }\n"
+                "}";
+
+    auto dbContext = DatabaseContext::create();
+    auto obj = createFromUserInput("EPSG:4326", dbContext);
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto cs = crs->coordinateSystem();
+    ASSERT_TRUE(cs != nullptr);
+    EXPECT_EQ(cs->exportToJSON(&(JSONFormatter::create()->setSchema("foo"))),
+              json);
+}
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 753878c5..8ee5814b 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -65,7 +65,8 @@ namespace {
 struct UnrelatedObject : public IComparable {
     UnrelatedObject() = default;
 
-    bool _isEquivalentTo(const IComparable *, Criterion) const override {
+    bool _isEquivalentTo(const IComparable *, Criterion,
+                         const DatabaseContextPtr &) const override {
         assert(false);
         return false;
     }
@@ -3257,7 +3258,9 @@ TEST(operation, webmerc_import_from_WKT2_EPSG_3785_deprecated) {
 
     EXPECT_EQ(
         crs->exportToWKT(
-            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL,
+                                 DatabaseContext::create())
+                .get()),
         "PROJCS[\"Popular Visualisation CRS / Mercator\",\n"
         "    GEOGCS[\"Popular Visualisation CRS\",\n"
         "        DATUM[\"Popular_Visualisation_Datum\",\n"
@@ -4628,6 +4631,33 @@ TEST(operation, geogCRS_to_geogCRS_context_WGS84_G1674_to_WGS84_G1762) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, geogCRS_to_geogCRS_context_EPSG_4240_Indian1975_to_EPSG_4326) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        authFactory->createCoordinateReferenceSystem("4240"), // Indian 1975
+        authFactory->createCoordinateReferenceSystem("4326"), ctxt);
+    ASSERT_EQ(list.size(), 3U);
+
+    // Indian 1975 to WGS 84 (4), 3.0 m, Thailand - onshore
+    EXPECT_EQ(list[0]->getEPSGCode(), 1812);
+
+    // The following is the one we want to see. It has a lesser accuracy than
+    // the above one and the same bbox, but the name of its area of use is
+    // slightly different
+    // Indian 1975 to WGS 84 (2), 5.0 m, Thailand - onshore and Gulf of Thailand
+    EXPECT_EQ(list[1]->getEPSGCode(), 1304);
+
+    // Indian 1975 to WGS 84 (3), 1.0 m, Thailand - Bongkot field
+    EXPECT_EQ(list[2]->getEPSGCode(), 1537);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, vertCRS_to_geogCRS_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");
@@ -4687,6 +4717,19 @@ TEST(operation, vertCRS_to_geogCRS_context) {
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
             "+proj=affine +s33=-0.304800609601219");
     }
+    {
+        auto ctxt =
+            CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            // NZVD2016 height
+            authFactory->createCoordinateReferenceSystem("7839"),
+            // NZGD2000
+            authFactory->createCoordinateReferenceSystem("4959"), ctxt);
+        ASSERT_EQ(list.size(), 1U);
+        EXPECT_EQ(
+            list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+            "+proj=vgridshift +grids=nzgeoid2016.gtx +multiplier=1");
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -5576,6 +5619,69 @@ TEST(operation, projCRS_no_id_to_geogCRS_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, projCRS_3D_to_geogCRS_3D_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto wkt = "PROJCRS[\"NAD83(HARN) / Oregon GIC Lambert (ft)\",\n"
+               "    BASEGEOGCRS[\"NAD83(HARN)\",\n"
+               "        DATUM[\"NAD83 (High Accuracy Reference Network)\",\n"
+               "            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
+               "                LENGTHUNIT[\"metre\",1]]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+               "        ID[\"EPSG\",4957]],\n"
+               "    CONVERSION[\"unnamed\",\n"
+               "        METHOD[\"Lambert Conic Conformal (2SP)\",\n"
+               "            ID[\"EPSG\",9802]],\n"
+               "        PARAMETER[\"Latitude of false origin\",41.75,\n"
+               "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+               "            ID[\"EPSG\",8821]],\n"
+               "        PARAMETER[\"Longitude of false origin\",-120.5,\n"
+               "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+               "            ID[\"EPSG\",8822]],\n"
+               "        PARAMETER[\"Latitude of 1st standard parallel\",43,\n"
+               "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+               "            ID[\"EPSG\",8823]],\n"
+               "        PARAMETER[\"Latitude of 2nd standard parallel\",45.5,\n"
+               "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+               "            ID[\"EPSG\",8824]],\n"
+               "        PARAMETER[\"Easting at false origin\",1312335.958,\n"
+               "            LENGTHUNIT[\"foot\",0.3048],\n"
+               "            ID[\"EPSG\",8826]],\n"
+               "        PARAMETER[\"Northing at false origin\",0,\n"
+               "            LENGTHUNIT[\"foot\",0.3048],\n"
+               "            ID[\"EPSG\",8827]]],\n"
+               "    CS[Cartesian,3],\n"
+               "        AXIS[\"easting\",east,\n"
+               "            ORDER[1],\n"
+               "            LENGTHUNIT[\"foot\",0.3048]],\n"
+               "        AXIS[\"northing\",north,\n"
+               "            ORDER[2],\n"
+               "            LENGTHUNIT[\"foot\",0.3048]],\n"
+               "        AXIS[\"ellipsoidal height (h)\",up,\n"
+               "            ORDER[3],\n"
+               "            LENGTHUNIT[\"foot\",0.3048]]]";
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto src = NN_CHECK_ASSERT(nn_dynamic_pointer_cast<CRS>(obj));
+    auto dst = authFactory->createCoordinateReferenceSystem(
+        "4957"); // NAD83(HARN) (3D)
+    auto list =
+        CoordinateOperationFactory::create()->createOperations(src, dst, ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              // Check that z ft->m conversion is done (and just once)
+              "+step +proj=unitconvert +xy_in=ft +z_in=ft +xy_out=m +z_out=m "
+              "+step +inv +proj=lcc +lat_0=41.75 +lon_0=-120.5 +lat_1=43 "
+              "+lat_2=45.5 +x_0=399999.9999984 +y_0=0 +ellps=GRS80 "
+              "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+              "+step +proj=axisswap +order=2,1");
+}
+// ---------------------------------------------------------------------------
+
 TEST(operation, projCRS_3D_to_projCRS_2D_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");
@@ -6509,6 +6615,18 @@ TEST(operation, transformation_VERTCON_to_PROJ_string) {
                   PROJStringFormatter::create().get()),
               "+proj=vgridshift +grids=bla.gtx +multiplier=0.001");
 }
+// ---------------------------------------------------------------------------
+
+TEST(operation, transformation_NZLVD_to_PROJ_string) {
+    auto dbContext = DatabaseContext::create();
+    auto factory = AuthorityFactory::create(dbContext, "EPSG");
+    auto op = factory->createCoordinateOperation("7860", false);
+    EXPECT_EQ(op->exportToPROJString(
+                  PROJStringFormatter::create(
+                      PROJStringFormatter::Convention::PROJ_5, dbContext)
+                      .get()),
+              "+proj=vgridshift +grids=auckht1946-nzvd2016.gtx +multiplier=1");
+}
 
 // ---------------------------------------------------------------------------
 
@@ -7134,6 +7252,22 @@ TEST(operation, vertCRS_to_vertCRS_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, vertCRS_to_vertCRS_New_Zealand_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // NZVD2016 height
+        authFactory->createCoordinateReferenceSystem("7839"),
+        // Auckland 1946 height
+        authFactory->createCoordinateReferenceSystem("5759"), ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=vgridshift +grids=auckht1946-nzvd2016.gtx +multiplier=1");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_to_geogCRS_3D) {
 
     auto compoundcrs_ft_obj = PROJStringParser().createFromPROJString(