Merge remote-tracking branch 'osgeo/master'

22 files changed, 1634 insertions, 1110 deletions

diff --git a/test/cli/Makefile.am b/test/cli/Makefile.am
index 4d11eaf3..46d9d36c 100644
--- a/test/cli/Makefile.am
+++ b/test/cli/Makefile.am
@@ -35,13 +35,13 @@ testprojinfo-check:
 	PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(TESTPROJINFO) $(PROJINFOEXE)
 
 test27-check:
-	$(TEST27) $(PROJEXE)
+	PROJ_LIB=$(PROJ_LIB) $(TEST27) $(PROJEXE)
 
 test83-check:
-	$(TEST83) $(PROJEXE)
+	PROJ_LIB=$(PROJ_LIB) $(TEST83) $(PROJEXE)
 
 testproj-check:
-	$(TESTPROJ) $(PROJEXE)
+	PROJ_LIB=$(PROJ_LIB) $(TESTPROJ) $(PROJEXE)
 
 testvarious-check:
 	PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(TESTVARIOUS) $(CS2CSEXE)
diff --git a/test/cli/td_out.dist b/test/cli/td_out.dist
index ab0c0911..6bad8e57 100644
--- a/test/cli/td_out.dist
+++ b/test/cli/td_out.dist
@@ -7,6 +7,18 @@ As above, but without ntv1 everything goes through conus file.
 111d00'00.000"W 44d00'00.000"N 0.0	111d0'2.788"W	43d59'59.725"N 0.000
 111d00'00.000"W 39d00'00.000"N 0.0	111d0'2.604"W	38d59'59.912"N 0.000
 ##############################################################
+Test --area Canada NAD27 NAD83 (using ntv1_can)
+43d59'59.732"N	111d0'3.208"W 0.000
+*	* inf
+##############################################################
+Test --bbox -141.01,40.04,-47.74,86.46 NAD27 NAD83 (using ntv1_can)
+43d59'59.732"N	111d0'3.208"W 0.000
+*	* inf
+##############################################################
+Test --area "USA - CONUS - onshore" NAD27 NAD83 (using conus)
+43d59'59.725"N	111d0'2.788"W 0.000
+38d59'59.912"N	111d0'2.604"W 0.000
+##############################################################
 Test MD used where available
 79d58'00.000"W 37d02'00.000"N 0.0	79d58'0.005"W	37d1'59.998"N 0.000
 79d58'00.000"W 36d58'00.000"N 0.0	79d57'59.128"W	36d58'0.501"N 0.000
diff --git a/test/cli/testcct b/test/cli/testcct
index 3fb0dd95..686931ea 100755
--- a/test/cli/testcct
+++ b/test/cli/testcct
@@ -32,6 +32,133 @@ echo "Testing cct -d 8 +proj=merc +R=1" >> ${OUT}
 echo "90 45" 0 | $EXE -d 8 +proj=merc +R=1 >>${OUT}
 echo "" >>${OUT}
 
+# tests without specifying the number of decimals (by default: 10 for radians and degrees, 4 for meters)
+echo "Testing echo 0.5 2 | cct -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad" >> ${OUT}
+echo 0.5 2 | $EXE -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad >> ${OUT}
+
+echo "Testing echo 0.5 2 | cct -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=deg" >> ${OUT}
+echo 0.5 2 | $EXE -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=deg >> ${OUT}
+
+echo "Testing echo 0.5 2 | cct -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=m +xy_out=km" >> ${OUT}
+echo 0.5 2 | $EXE -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=m +xy_out=km >> ${OUT}
+echo "" >> ${OUT}
+
+# tests for which the number of decimals has been specified (-d 6)
+echo "Testing echo 0.5 2 | cct -d 6 -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad" >> ${OUT}
+echo 0.5 2 | $EXE -d 6 -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad >> ${OUT}
+
+echo "Testing echo 0.5 2 | cct -d 6 -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=deg" >> ${OUT}
+echo 0.5 2 | $EXE -d 6 -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=deg >> ${OUT}
+
+echo "Testing echo 0.5 2 | cct -d 6 -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=m +xy_out=km" >> ${OUT}
+echo 0.5 2 | $EXE -d 6 -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=m +xy_out=km >> ${OUT}
+echo "" >> ${OUT}
+
+echo "Test cct with object code initialization" >> ${OUT}
+echo "3541657.3778 948984.2343 5201383.5231 2020.5" | $EXE EPSG:8366 >>${OUT}
+
+echo "Test cct with object name initialization" >> ${OUT}
+echo "3541657.3778 948984.2343 5201383.5231 2020.5" | $EXE "ITRF2014 to ETRF2014 (1)"  >>${OUT}
+
+echo "Test cct with object code initialization and file input" >> ${OUT}
+echo "3541657.3778 948984.2343 5201383.5231 2020.5" > a
+echo "3541658.0000 948985.0000 5201384.0000 2020.5" > b
+$EXE EPSG:8366 a b >>${OUT}
+/bin/rm a b
+
+cat > in.wkt <<EOF
+COORDINATEOPERATION["ITRF2014 to ETRF2014 (1)",
+    VERSION["EUREF-Eur"],
+    SOURCECRS[
+        GEODCRS["ITRF2014",
+            DYNAMIC[
+                FRAMEEPOCH[2010]],
+            DATUM["International Terrestrial Reference Frame 2014",
+                ELLIPSOID["GRS 1980",6378137,298.257222101,
+                    LENGTHUNIT["metre",1]]],
+            PRIMEM["Greenwich",0,
+                ANGLEUNIT["degree",0.0174532925199433]],
+            CS[Cartesian,3],
+                AXIS["(X)",geocentricX,
+                    ORDER[1],
+                    LENGTHUNIT["metre",1]],
+                AXIS["(Y)",geocentricY,
+                    ORDER[2],
+                    LENGTHUNIT["metre",1]],
+                AXIS["(Z)",geocentricZ,
+                    ORDER[3],
+                    LENGTHUNIT["metre",1]],
+            ID["EPSG",7789]]],
+    TARGETCRS[
+        GEODCRS["ETRF2014",
+            DATUM["European Terrestrial Reference Frame 2014",
+                ELLIPSOID["GRS 1980",6378137,298.257222101,
+                    LENGTHUNIT["metre",1]]],
+            PRIMEM["Greenwich",0,
+                ANGLEUNIT["degree",0.0174532925199433]],
+            CS[Cartesian,3],
+                AXIS["(X)",geocentricX,
+                    ORDER[1],
+                    LENGTHUNIT["metre",1]],
+                AXIS["(Y)",geocentricY,
+                    ORDER[2],
+                    LENGTHUNIT["metre",1]],
+                AXIS["(Z)",geocentricZ,
+                    ORDER[3],
+                    LENGTHUNIT["metre",1]],
+            ID["EPSG",8401]]],
+    METHOD["Time-dependent Position Vector tfm (geocentric)",
+        ID["EPSG",1053]],
+    PARAMETER["X-axis translation",0,
+        LENGTHUNIT["millimetre",0.001],
+        ID["EPSG",8605]],
+    PARAMETER["Y-axis translation",0,
+        LENGTHUNIT["millimetre",0.001],
+        ID["EPSG",8606]],
+    PARAMETER["Z-axis translation",0,
+        LENGTHUNIT["millimetre",0.001],
+        ID["EPSG",8607]],
+    PARAMETER["X-axis rotation",0,
+        ANGLEUNIT["milliarc-second",4.84813681109536E-09],
+        ID["EPSG",8608]],
+    PARAMETER["Y-axis rotation",0,
+        ANGLEUNIT["milliarc-second",4.84813681109536E-09],
+        ID["EPSG",8609]],
+    PARAMETER["Z-axis rotation",0,
+        ANGLEUNIT["milliarc-second",4.84813681109536E-09],
+        ID["EPSG",8610]],
+    PARAMETER["Scale difference",0,
+        SCALEUNIT["parts per billion",1E-09],
+        ID["EPSG",8611]],
+    PARAMETER["Rate of change of X-axis translation",0,
+        LENGTHUNIT["millimetres per year",3.16887651727315E-11],
+        ID["EPSG",1040]],
+    PARAMETER["Rate of change of Y-axis translation",0,
+        LENGTHUNIT["millimetres per year",3.16887651727315E-11],
+        ID["EPSG",1041]],
+    PARAMETER["Rate of change of Z-axis translation",0,
+        LENGTHUNIT["millimetres per year",3.16887651727315E-11],
+        ID["EPSG",1042]],
+    PARAMETER["Rate of change of X-axis rotation",0.085,
+        ANGLEUNIT["milliarc-seconds per year",1.53631468932076E-16],
+        ID["EPSG",1043]],
+    PARAMETER["Rate of change of Y-axis rotation",0.531,
+        ANGLEUNIT["milliarc-seconds per year",1.53631468932076E-16],
+        ID["EPSG",1044]],
+    PARAMETER["Rate of change of Z-axis rotation",-0.77,
+        ANGLEUNIT["milliarc-seconds per year",1.53631468932076E-16],
+        ID["EPSG",1045]],
+    PARAMETER["Rate of change of Scale difference",0,
+        SCALEUNIT["parts per billion per year",3.16887651727315E-17],
+        ID["EPSG",1046]],
+    PARAMETER["Parameter reference epoch",1989,
+        TIMEUNIT["year",31556925.445],
+        ID["EPSG",1047]]]
+EOF
+echo "Test cct with WKT in a file" >> ${OUT}
+echo "3541657.3778 948984.2343 5201383.5231 2020.5" | $EXE @in.wkt >>${OUT}
+rm in.wkt
+
 # do 'diff' with distribution results
 echo "diff ${OUT} with testcct_out.dist"
 diff -u ${OUT} ${TEST_CLI_DIR}/testcct_out.dist
diff --git a/test/cli/testcct_out.dist b/test/cli/testcct_out.dist
index 44dd6964..7762ace7 100644
--- a/test/cli/testcct_out.dist
+++ b/test/cli/testcct_out.dist
@@ -1,3 +1,26 @@
 Testing cct -d 8 +proj=merc +R=1
    1.57079633     0.88137359    0.00000000           inf
 
+Testing echo 0.5 2 | cct -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad
+  0.5000000000    2.0000000000        0.0000        0.0000
+Testing echo 0.5 2 | cct -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=deg
+  0.5000000000    2.0000000000        0.0000        0.0000
+Testing echo 0.5 2 | cct -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=m +xy_out=km
+       0.0005         0.0020        0.0000        0.0000
+
+Testing echo 0.5 2 | cct -d 6 -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad
+      0.500000        2.000000      0.000000        0.0000
+Testing echo 0.5 2 | cct -d 6 -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=deg
+      0.500000        2.000000      0.000000        0.0000
+Testing echo 0.5 2 | cct -d 6 -z 0 -t 0 +proj=pipeline +step +proj=unitconvert +xy_in=m +xy_out=km
+     0.000500       0.002000      0.000000        0.0000
+
+Test cct with object code initialization
+ 3541657.9112    948983.7503  5201383.2482     2020.5000
+Test cct with object name initialization
+ 3541657.9112    948983.7503  5201383.2482     2020.5000
+Test cct with object code initialization and file input
+ 3541657.9112    948983.7503  5201383.2482     2020.5000
+ 3541658.5334    948984.5160  5201383.7251     2020.5000
+Test cct with WKT in a file
+ 3541657.9112    948983.7503  5201383.2482     2020.5000
diff --git a/test/cli/testdatumfile b/test/cli/testdatumfile
index 16e4bbc3..5a013f12 100755
--- a/test/cli/testdatumfile
+++ b/test/cli/testdatumfile
@@ -60,6 +60,27 @@ $EXE +proj=latlong +ellps=clrk66 '+nadgrids="./dir with "" space/myconus"' \
 EOF
 
 echo "##############################################################" >> ${OUT}
+echo "Test --area Canada NAD27 NAD83 (using ntv1_can)" >> ${OUT}
+$EXE --area Canada NAD27 NAD83 >>${OUT} <<EOF
+44 -111
+39 -111
+EOF
+
+echo "##############################################################" >> ${OUT}
+echo "Test --bbox -141.01,40.04,-47.74,86.46 NAD27 NAD83 (using ntv1_can)" >> ${OUT}
+$EXE --bbox -141.01,40.04,-47.74,86.46 NAD27 NAD83 >>${OUT} <<EOF
+44 -111
+39 -111
+EOF
+
+echo "##############################################################" >> ${OUT}
+echo "Test --area \"USA - CONUS - onshore\" NAD27 NAD83 (using conus)" >> ${OUT}
+$EXE --area "USA - CONUS - onshore" NAD27 NAD83 >>${OUT} <<EOF
+44 -111
+39 -111
+EOF
+
+echo "##############################################################" >> ${OUT}
 echo Test MD used where available >> ${OUT}
 #
 $EXE +proj=latlong +ellps=clrk66 +nadgrids=MD,conus \
diff --git a/test/cli/testprojinfo b/test/cli/testprojinfo
index 62713af5..c31cfef0 100755
--- a/test/cli/testprojinfo
+++ b/test/cli/testprojinfo
@@ -171,6 +171,10 @@ echo 'Testing -k datum EPSG:6326' >> ${OUT}
 $EXE -k datum EPSG:6326 >>${OUT} 2>&1
 echo "" >>${OUT}
 
+echo 'Testing -k ensemble WGS84' >> ${OUT}
+$EXE -k ensemble WGS84 >>${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}
diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
index 245bb258..8e8ef294 100644
--- a/test/cli/testprojinfo_out.dist
+++ b/test/cli/testprojinfo_out.dist
@@ -4,9 +4,16 @@ PROJ.4 string:
 
 WKT2:2019 string:
 GEOGCRS["WGS 84",
-    DATUM["World Geodetic System 1984",
+    ENSEMBLE["World Geodetic System 1984 ensemble",
+        MEMBER["World Geodetic System 1984 (Transit)"],
+        MEMBER["World Geodetic System 1984 (G730)"],
+        MEMBER["World Geodetic System 1984 (G873)"],
+        MEMBER["World Geodetic System 1984 (G1150)"],
+        MEMBER["World Geodetic System 1984 (G1674)"],
+        MEMBER["World Geodetic System 1984 (G1762)"],
         ELLIPSOID["WGS 84",6378137,298.257223563,
-            LENGTHUNIT["metre",1]]],
+            LENGTHUNIT["metre",1]],
+        ENSEMBLEACCURACY[2.0]],
     PRIMEM["Greenwich",0,
         ANGLEUNIT["degree",0.0174532925199433]],
     CS[ellipsoidal,2],
@@ -59,9 +66,16 @@ GEODCRS["WGS 84",
 Testing projinfo -o WKT2_2019 EPSG:4326
 WKT2:2019 string:
 GEOGCRS["WGS 84",
-    DATUM["World Geodetic System 1984",
+    ENSEMBLE["World Geodetic System 1984 ensemble",
+        MEMBER["World Geodetic System 1984 (Transit)"],
+        MEMBER["World Geodetic System 1984 (G730)"],
+        MEMBER["World Geodetic System 1984 (G873)"],
+        MEMBER["World Geodetic System 1984 (G1150)"],
+        MEMBER["World Geodetic System 1984 (G1674)"],
+        MEMBER["World Geodetic System 1984 (G1762)"],
         ELLIPSOID["WGS 84",6378137,298.257223563,
-            LENGTHUNIT["metre",1]]],
+            LENGTHUNIT["metre",1]],
+        ENSEMBLEACCURACY[2.0]],
     PRIMEM["Greenwich",0,
         ANGLEUNIT["degree",0.0174532925199433]],
     CS[ellipsoidal,2],
@@ -102,9 +116,16 @@ GEODCRS["WGS 84",
 
 WKT2:2019 string:
 GEOGCRS["WGS 84",
-    DATUM["World Geodetic System 1984",
+    ENSEMBLE["World Geodetic System 1984 ensemble",
+        MEMBER["World Geodetic System 1984 (Transit)"],
+        MEMBER["World Geodetic System 1984 (G730)"],
+        MEMBER["World Geodetic System 1984 (G873)"],
+        MEMBER["World Geodetic System 1984 (G1150)"],
+        MEMBER["World Geodetic System 1984 (G1674)"],
+        MEMBER["World Geodetic System 1984 (G1762)"],
         ELLIPSOID["WGS 84",6378137,298.257223563,
-            LENGTHUNIT["metre",1]]],
+            LENGTHUNIT["metre",1]],
+        ENSEMBLEACCURACY[2.0]],
     PRIMEM["Greenwich",0,
         ANGLEUNIT["degree",0.0174532925199433]],
     CS[ellipsoidal,2],
@@ -142,13 +163,61 @@ PROJJSON:
   "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json",
   "type": "GeographicCRS",
   "name": "WGS 84",
-  "datum": {
-    "type": "GeodeticReferenceFrame",
-    "name": "World Geodetic System 1984",
+  "datum_ensemble": {
+    "name": "World Geodetic System 1984 ensemble",
+    "members": [
+      {
+        "name": "World Geodetic System 1984 (Transit)",
+        "id": {
+          "authority": "EPSG",
+          "code": 1166
+        }
+      },
+      {
+        "name": "World Geodetic System 1984 (G730)",
+        "id": {
+          "authority": "EPSG",
+          "code": 1152
+        }
+      },
+      {
+        "name": "World Geodetic System 1984 (G873)",
+        "id": {
+          "authority": "EPSG",
+          "code": 1153
+        }
+      },
+      {
+        "name": "World Geodetic System 1984 (G1150)",
+        "id": {
+          "authority": "EPSG",
+          "code": 1154
+        }
+      },
+      {
+        "name": "World Geodetic System 1984 (G1674)",
+        "id": {
+          "authority": "EPSG",
+          "code": 1155
+        }
+      },
+      {
+        "name": "World Geodetic System 1984 (G1762)",
+        "id": {
+          "authority": "EPSG",
+          "code": 1156
+        }
+      }
+    ],
     "ellipsoid": {
       "name": "WGS 84",
       "semi_major_axis": 6378137,
       "inverse_flattening": 298.257223563
+    },
+    "accuracy": "2.0",
+    "id": {
+      "authority": "EPSG",
+      "code": 6326
     }
   },
   "coordinate_system": {
@@ -317,7 +386,8 @@ COORDINATEOPERATION["NAD27 to NAD83 (3)",
         SCOPE["Historic record only - now superseded - see remarks."],
         AREA["Canada - onshore and offshore - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon."],
         BBOX[40.04,-141.01,86.46,-47.74]],
-    ID["DERIVED_FROM(EPSG)",1312]]
+    ID["DERIVED_FROM(EPSG)",1312],
+    REMARK["Uses NTv1 method. Replaced in Quebec by code 1462 and elsewhere in 1997 by NTv2 (transformation code 1313). Input expects longitudes to be positive west; EPSG GeogCRS NAD27 (code 4267) and NAD83 (code 4269) have longitudes positive east."]]
 
 -------------------------------------
 Operation No. 2:
@@ -371,7 +441,8 @@ COORDINATEOPERATION["NAD27 to NAD83 (4)",
         SCOPE["Transformation of coordinates at 1m to 2m level of accuracy."],
         AREA["Canada - onshore - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon; offshore east coast."],
         BBOX[40.04,-141.01,83.17,-47.74]],
-    ID["DERIVED_FROM(EPSG)",1313]]
+    ID["DERIVED_FROM(EPSG)",1313],
+    REMARK["Uses NTv2 data files. Replaces NTv1 (transformation code 1312) except in Quebec. Input expects longitudes to be positive west; EPSG GeogCRS NAD27 (code 4267) and (code 4269) have longitudes positive east. May be used as tfm to WGS 84 - see code 1693."]]
 
 -------------------------------------
 Operation No. 3:
@@ -425,7 +496,8 @@ COORDINATEOPERATION["NAD27 to NAD83 (1)",
         SCOPE["Transformation of coordinates at 0.2m level of accuracy."],
         AREA["United States (USA) - CONUS including EEZ -onshore and offshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Gulf of Mexico (GoM) OCS."],
         BBOX[23.81,-129.17,49.38,-65.69]],
-    ID["DERIVED_FROM(EPSG)",1241]]
+    ID["DERIVED_FROM(EPSG)",1241],
+    REMARK["Uses NADCON method which expects longitudes positive west; EPSG GeogCRS NAD27 (code 4267) and NAD83 (code 4269) have longitudes positive east."]]
 
 -------------------------------------
 Operation No. 4:
@@ -479,7 +551,8 @@ COORDINATEOPERATION["NAD27 to NAD83 (2)",
         SCOPE["Geodesy."],
         AREA["United States (USA) - Alaska including EEZ."],
         BBOX[47.88,167.65,74.71,-129.99]],
-    ID["DERIVED_FROM(EPSG)",1243]]
+    ID["DERIVED_FROM(EPSG)",1243],
+    REMARK["Uses NADCON method which expects longitudes positive west; EPSG GeogCRS NAD27 (code 4267) and NAD83 (code 4269) have longitudes positive east. May be used as transformation to WGS 84 - see NAD27 to WGS 84 (85) (code 15864)."]]
 
 -------------------------------------
 Operation No. 5:
@@ -533,7 +606,8 @@ COORDINATEOPERATION["NAD27 to NAD83 (6)",
         SCOPE["Transformation of coordinates at 1m to 2m level of accuracy."],
         AREA["Canada - Quebec."],
         BBOX[44.99,-79.85,62.62,-57.1]],
-    ID["DERIVED_FROM(EPSG)",1573]]
+    ID["DERIVED_FROM(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."]]
 
 -------------------------------------
 Operation No. 6:
@@ -916,7 +990,7 @@ PROJ.4 string:
 +proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs
 
 Testing RH2000 height to SWEREF99: projinfo -s EPSG:5613 -t EPSG:4977
-Candidate operations found: 1
+Candidate operations found: 2
 -------------------------------------
 Operation No. 1:
 
@@ -968,8 +1042,55 @@ COORDINATEOPERATION["Inverse of SWEREF99 to RH2000 height",
         BBOX[55.28,10.93,69.07,24.17]],
     ID["INVERSE(DERIVED_FROM(PROJ))","EPSG_4977_TO_EPSG_5613"]]
 
+-------------------------------------
+Operation No. 2:
+
+unknown id, Transformation from RH2000 height to SWEREF99 (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation
+
+PROJ string:
++proj=noop
+
+WKT2:2019 string:
+COORDINATEOPERATION["Transformation from RH2000 height to SWEREF99 (ballpark vertical transformation, without ellipsoid height to vertical height correction)",
+    SOURCECRS[
+        VERTCRS["RH2000 height",
+            DYNAMIC[
+                FRAMEEPOCH[2000]],
+            VDATUM["Rikets hojdsystem 2000"],
+            CS[vertical,1],
+                AXIS["gravity-related height (H)",up,
+                    LENGTHUNIT["metre",1]],
+            ID["EPSG",5613]]],
+    TARGETCRS[
+        GEOGCRS["SWEREF99",
+            DATUM["SWEREF99",
+                ELLIPSOID["GRS 1980",6378137,298.257222101,
+                    LENGTHUNIT["metre",1]]],
+            PRIMEM["Greenwich",0,
+                ANGLEUNIT["degree",0.0174532925199433]],
+            CS[ellipsoidal,3],
+                AXIS["geodetic latitude (Lat)",north,
+                    ORDER[1],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+                AXIS["geodetic longitude (Lon)",east,
+                    ORDER[2],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+                AXIS["ellipsoidal height (h)",up,
+                    ORDER[3],
+                    LENGTHUNIT["metre",1]],
+            ID["EPSG",4977]]],
+    METHOD["Change of Vertical Unit",
+        ID["EPSG",1069]],
+    PARAMETER["Unit conversion scalar",1,
+        SCALEUNIT["unity",1],
+        ID["EPSG",1051]],
+    USAGE[
+        SCOPE["unknown"],
+        AREA["World"],
+        BBOX[-90,-180,90,180]]]
+
 Testing NAD83(2011) + NAVD88 height -> NAD83(2011) : projinfo -s EPSG:6349 -t EPSG:6319 --spatial-test intersects -o PROJ
-Candidate operations found: 2
+Candidate operations found: 3
 -------------------------------------
 Operation No. 1:
 
@@ -996,8 +1117,16 @@ PROJ string:
   +step +proj=unitconvert +xy_in=rad +xy_out=deg
   +step +proj=axisswap +order=2,1
 
+-------------------------------------
+Operation No. 3:
+
+unknown id, Transformation from NAVD88 height to NAD83(2011) (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation
+
+PROJ string:
++proj=noop
+
 Testing NGF IGN69 height to RGF93: projinfo -s EPSG:5720 -t EPSG:4965 -o PROJ
-Candidate operations found: 1
+Candidate operations found: 2
 -------------------------------------
 Operation No. 1:
 
@@ -1011,6 +1140,14 @@ PROJ string:
   +step +proj=unitconvert +xy_in=rad +xy_out=deg
   +step +proj=axisswap +order=2,1
 
+-------------------------------------
+Operation No. 2:
+
+unknown id, Transformation from NGF-IGN69 height to RGF93 (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation
+
+PROJ string:
++proj=noop
+
 Testing EPSG:32631 --3d
 PROJ.4 string:
 +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +type=crs
@@ -1018,9 +1155,16 @@ PROJ.4 string:
 WKT2:2019 string:
 PROJCRS["WGS 84 / UTM zone 31N",
     BASEGEOGCRS["WGS 84",
-        DATUM["World Geodetic System 1984",
+        ENSEMBLE["World Geodetic System 1984 ensemble",
+            MEMBER["World Geodetic System 1984 (Transit)"],
+            MEMBER["World Geodetic System 1984 (G730)"],
+            MEMBER["World Geodetic System 1984 (G873)"],
+            MEMBER["World Geodetic System 1984 (G1150)"],
+            MEMBER["World Geodetic System 1984 (G1674)"],
+            MEMBER["World Geodetic System 1984 (G1762)"],
             ELLIPSOID["WGS 84",6378137,298.257223563,
-                LENGTHUNIT["metre",1]]],
+                LENGTHUNIT["metre",1]],
+            ENSEMBLEACCURACY[2.0]],
         PRIMEM["Greenwich",0,
             ANGLEUNIT["degree",0.0174532925199433]],
         ID["EPSG",4979]],
@@ -1179,6 +1323,27 @@ DATUM["World Geodetic System 1984",
         LENGTHUNIT["metre",1]],
     ID["EPSG",6326]]
 
+Testing -k ensemble WGS84
+WKT2:2019 string:
+ENSEMBLE["World Geodetic System 1984 ensemble",
+    MEMBER["World Geodetic System 1984 (Transit)",
+        ID["EPSG",1166]],
+    MEMBER["World Geodetic System 1984 (G730)",
+        ID["EPSG",1152]],
+    MEMBER["World Geodetic System 1984 (G873)",
+        ID["EPSG",1153]],
+    MEMBER["World Geodetic System 1984 (G1150)",
+        ID["EPSG",1154]],
+    MEMBER["World Geodetic System 1984 (G1674)",
+        ID["EPSG",1155]],
+    MEMBER["World Geodetic System 1984 (G1762)",
+        ID["EPSG",1156]],
+    ELLIPSOID["WGS 84",6378137,298.257223563,
+        LENGTHUNIT["metre",1],
+        ID["EPSG",7030]],
+    ENSEMBLEACCURACY[2.0],
+    ID["EPSG",6326]]
+
 Testing -k operation EPSG:8457 -o PROJ -q
 +proj=pipeline
   +step +proj=axisswap +order=2,1
@@ -1193,9 +1358,23 @@ Testing -k operation EPSG:8457 -o PROJ -q
 
 Testing D_WGS_1984
 WKT2:2019 string:
-DATUM["World Geodetic System 1984",
+ENSEMBLE["World Geodetic System 1984 ensemble",
+    MEMBER["World Geodetic System 1984 (Transit)",
+        ID["EPSG",1166]],
+    MEMBER["World Geodetic System 1984 (G730)",
+        ID["EPSG",1152]],
+    MEMBER["World Geodetic System 1984 (G873)",
+        ID["EPSG",1153]],
+    MEMBER["World Geodetic System 1984 (G1150)",
+        ID["EPSG",1154]],
+    MEMBER["World Geodetic System 1984 (G1674)",
+        ID["EPSG",1155]],
+    MEMBER["World Geodetic System 1984 (G1762)",
+        ID["EPSG",1156]],
     ELLIPSOID["WGS 84",6378137,298.257223563,
-        LENGTHUNIT["metre",1]],
+        LENGTHUNIT["metre",1],
+        ID["EPSG",7030]],
+    ENSEMBLEACCURACY[2.0],
     ID["EPSG",6326]]
 
 Testing -k datum D_WGS_1984
diff --git a/test/cli/testvarious b/test/cli/testvarious
index 292ee316..82be4992 100755
--- a/test/cli/testvarious
+++ b/test/cli/testvarious
@@ -1009,6 +1009,19 @@ $EXE -f %.3f EPSG:4686 EPSG:6247 -E >> ${OUT} <<EOF
 4.8 -74.25
 EOF
 
+echo  "##############################################################" >> ${OUT}
+echo  "Test effect of --authority (https://github.com/OSGeo/PROJ/issues/2442)" >> ${OUT}
+echo  "This test might be a bit fragile if proj.db content changes" >> ${OUT}
+echo  "The first result should use the 'WGS_1984_(ITRF08)_To_NAD_1983_2011' (ESRI:108363) operation" >> ${OUT}
+echo  "and the second one a no-op" >> ${OUT}
+$EXE -E +proj=latlong +datum=WGS84 +no_defs +to +init=epsg:6342 >> ${OUT} <<EOF
+-105 40
+EOF
+$EXE --authority EPSG -E +proj=latlong +datum=WGS84 +no_defs +to +init=epsg:6342 >> ${OUT} <<EOF
+-105 40
+EOF
+
+
 # Done!
 # do 'diff' with distribution results
 echo "diff ${OUT} with ${OUT}.dist"
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index 70b2ab6e..fe1aa452 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -485,3 +485,10 @@ Test EPSG:xxxx EPSG:yyyy filename
 ##############################################################
 Test Colombia Urban
 4.8 -74.25	122543.174	80859.033 0.000
+##############################################################
+Test effect of --authority (https://github.com/OSGeo/PROJ/issues/2442)
+This test might be a bit fragile if proj.db content changes
+The first result should use the 'WGS_1984_(ITRF08)_To_NAD_1983_2011' (ESRI:108363) operation
+and the second one a no-op
+-105 40	500000.86	4427756.50 0.00
+-105 40	500000.00	4427757.22 0.00
diff --git a/test/fuzzers/build.sh b/test/fuzzers/build.sh
index 1b2bf79c..37a207e0 100755
--- a/test/fuzzers/build.sh
+++ b/test/fuzzers/build.sh
@@ -69,12 +69,8 @@ build_fuzzer()
         $LIB_FUZZING_ENGINE src/.libs/libproj.a $EXTRA_LIBS
 }
 
-build_fuzzer standard_fuzzer test/fuzzers/standard_fuzzer.cpp
 build_fuzzer proj_crs_to_crs_fuzzer test/fuzzers/proj_crs_to_crs_fuzzer.cpp
 
-echo "[libfuzzer]" > $OUT/standard_fuzzer.options
-echo "max_len = 10000" >> $OUT/standard_fuzzer.options
-
 echo "[libfuzzer]" > $OUT/proj_crs_to_crs_fuzzer.options
 echo "max_len = 10000" >> $OUT/proj_crs_to_crs_fuzzer.options
 
diff --git a/test/fuzzers/standard_fuzzer.cpp b/test/fuzzers/standard_fuzzer.cpp
deleted file mode 100644
index 468e8cbb..00000000
--- a/test/fuzzers/standard_fuzzer.cpp
+++ /dev/null
@@ -1,218 +0,0 @@
-/******************************************************************************
- *
- * Project:  proj.4
- * Purpose:  Fuzzer
- * Author:   Even Rouault, even.rouault at spatialys.com
- *
- ******************************************************************************
- * Copyright (c) 2017, Even Rouault <even.rouault at spatialys.com>
- *
- * 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.
- ****************************************************************************/
-
-#include <stddef.h>
-#include <stdint.h>
-#include <stdlib.h>
-#include <string.h>
-#include <stdio.h>
-#include <sys/types.h>
-#include <unistd.h>
-
-#define ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
-#include "proj.h"
-#include "proj_api.h"
-
-/* Standalone build:
-g++ -g -std=c++11 standard_fuzzer.cpp -o standard_fuzzer -fvisibility=hidden -DSTANDALONE ../../src/.libs/libproj.a -lpthread -lsqlite3 -I../../src -I../../include
-*/
-
-extern "C" int LLVMFuzzerInitialize(int* argc, char*** argv);
-extern "C" int LLVMFuzzerTestOneInput(const uint8_t *buf, size_t len);
-
-int LLVMFuzzerInitialize(int* /*argc*/, char*** argv)
-{
-    const char* argv0 = (*argv)[0];
-    char* path = pj_strdup(argv0);
-    char* lastslash = strrchr(path, '/');
-    if( lastslash )
-    {
-        *lastslash = '\0';
-        setenv("PROJ_LIB", path, 1);
-    }
-    else
-    {
-        setenv("PROJ_LIB", ".", 1);
-    }
-    free(path);
-    return 0;
-}
-
-int LLVMFuzzerTestOneInput(const uint8_t *buf, size_t len)
-{
-    if( len > 1000 )
-    {
-#ifdef STANDALONE
-        fprintf(stderr, "Input too large\n");
-#endif
-        return 0;
-    }
-
-    /* We expect the blob to be 3 lines: */
-    /* source proj string\ndestination proj string\nx y */
-    char* buf_dup = (char*)malloc(len+1);
-    memcpy(buf_dup, buf, len);
-    buf_dup[len] = 0;
-    char* first_line = buf_dup;
-    char* first_newline = strchr(first_line, '\n');
-    if( !first_newline )
-    {
-        free(buf_dup);
-        return 0;
-    }
-    first_newline[0] = 0;
-    char* second_line = first_newline + 1;
-    char* second_newline = strchr(second_line, '\n');
-    if( !second_newline )
-    {
-        free(buf_dup);
-        return 0;
-    }
-    second_newline[0] = 0;
-    char* third_line = second_newline + 1;
-    projPJ pj_src = pj_init_plus(first_line);
-    if( !pj_src )
-    {
-        free(buf_dup);
-        return 0;
-    }
-    projPJ pj_dst = pj_init_plus(second_line);
-    if( !pj_dst )
-    {
-        free(buf_dup);
-        pj_free(pj_src);
-        proj_cleanup();
-        return 0;
-    }
-    double x = 0, y = 0, z = 9;
-    bool from_binary = false;
-    bool has_z = false;
-    if( strncmp(third_line, "BINARY_2D:", strlen("BINARY_2D:")) == 0 &&
-        third_line - first_line + strlen("BINARY_2D:") + 2 * sizeof(double) <= len )
-    {
-        from_binary = true;
-        memcpy(&x, third_line + strlen("BINARY_2D:"), sizeof(double));
-        memcpy(&y, third_line + strlen("BINARY_2D:") + sizeof(double), sizeof(double));
-    }
-    else if( strncmp(third_line, "BINARY_3D:", strlen("BINARY_3D:")) == 0 &&
-             third_line - first_line + strlen("BINARY_3D:") + 3 * sizeof(double) <= len )
-    {
-        from_binary = true;
-        has_z = true;
-        memcpy(&x, third_line + strlen("BINARY_3D:"), sizeof(double));
-        memcpy(&y, third_line + strlen("BINARY_3D:") + sizeof(double), sizeof(double));
-        memcpy(&z, third_line + strlen("BINARY_3D:") + 2 * sizeof(double), sizeof(double));
-    }
-    else if( sscanf(third_line, "%lf %lf", &x, &y) != 2 )
-    {
-        free(buf_dup);
-        pj_free(pj_src);
-        pj_free(pj_dst);
-        proj_cleanup();
-        return 0;
-    }
-#ifdef STANDALONE
-    fprintf(stderr, "src=%s\n", first_line);
-    fprintf(stderr, "dst=%s\n", second_line);
-    if( from_binary )
-    {
-        if( has_z )
-            fprintf(stderr, "coord (from binary)=%.18g %.18g %.18g\n", x, y, z);
-        else
-            fprintf(stderr, "coord (from binary)=%.18g %.18g\n", x, y);
-    }
-    else
-        fprintf(stderr, "coord=%s\n", third_line);
-#endif
-    if( has_z )
-        pj_transform( pj_src, pj_dst, 1, 0, &x, &y, &z );
-    else
-        pj_transform( pj_src, pj_dst, 1, 0, &x, &y, NULL );
-    free(buf_dup);
-    pj_free(pj_src);
-    pj_free(pj_dst);
-    proj_cleanup();
-    return 0;
-}
-
-#ifdef STANDALONE
-
-int main(int argc, char* argv[])
-{
-    if( argc < 2 )
-    {
-        const char str[] =
-            "+proj=longlat +datum=WGS84 +nodefs\n+proj=longlat +datum=WGS84 +nodefs\n2 49";
-        int ret = LLVMFuzzerTestOneInput((const uint8_t*)(str), sizeof(str) - 1);
-        if( ret )
-            return ret;
-
-        const char str2[] =
-            "+proj=longlat +datum=WGS84 +nodefs\n+proj=longlat +datum=WGS84 +nodefs\nBINARY_2D:\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\xFF\xFF";
-        ret = LLVMFuzzerTestOneInput((const uint8_t*)(str2), sizeof(str2) - 1);
-        if( ret )
-            return ret;
-
-        const char str3[] =
-            "+proj=longlat +datum=WGS84 +nodefs\n+proj=longlat +datum=WGS84 +nodefs\nBINARY_3D:\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\xFF\xFF\x00\x00\x00\x00\x00\x00\x00\x00";
-        ret = LLVMFuzzerTestOneInput((const uint8_t*)(str3), sizeof(str3) - 1);
-        if( ret )
-            return ret;
-
-        return 0;
-    }
-    else
-    {
-        int nRet = 0;
-        void* buf = NULL;
-        int nLen = 0;
-        FILE* f = fopen(argv[1], "rb");
-        if( !f )
-        {
-            fprintf(stderr, "%s does not exist.\n", argv[1]);
-            exit(1);
-        }
-        fseek(f, 0, SEEK_END);
-        nLen = (int)ftell(f);
-        fseek(f, 0, SEEK_SET);
-        buf = malloc(nLen);
-        if( !buf )
-        {
-            fprintf(stderr, "malloc failed.\n");
-            fclose(f);
-            exit(1);
-        }
-        fread(buf, nLen, 1, f);
-        fclose(f);
-        nRet = LLVMFuzzerTestOneInput((const uint8_t*)(buf), nLen);
-        free(buf);
-        return nRet;
-    }
-}
-
-#endif // STANDALONE
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index cfce5041..813b67bb 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -1423,9 +1423,9 @@ expect  -0.001790143 -0.000895247
 ===============================================================================
 
 -------------------------------------------------------------------------------
-operation +proj=etmerc   +ellps=GRS80 +zone=30
+operation +proj=etmerc   +ellps=GRS80
 -------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 50 nm
 accept  2 1
 expect  222650.796797586 110642.229411933
 accept  2 -1
@@ -1434,17 +1434,28 @@ accept  -2 1
 expect  -222650.796797586 110642.229411933
 accept  -2 -1
 expect  -222650.796797586 -110642.229411933
+# near pole
+accept  30 89.9999
+expect  5.584698978 10001956.056248082
+# 3900 km from central meridian
+accept  44.69 35.37
+expect  4168136.489446198 4985511.302287407
 
 direction inverse
 accept  200 100
-expect  0.001796631 0.000904369
+expect  0.00179663056816 0.00090436947663
 accept  200 -100
-expect  0.001796631 -0.000904369
+expect  0.00179663056816 -0.00090436947663
 accept  -200 100
-expect  -0.001796631 0.000904369
+expect  -0.00179663056816 0.00090436947663
 accept  -200 -100
-expect  -0.001796631 -0.000904369
-
+expect  -0.00179663056816 -0.00090436947663
+# near pole
+accept  6 1.0001e7
+expect  0.35596960759234 89.99135362646302
+# 3900 km from central meridian
+accept  4168136.489446198 4985511.302287407
+expect  44.69 35.37
 
 ===============================================================================
 # Fahey
@@ -3355,30 +3366,54 @@ expect  -0.001953415 -0.000820580
 -------------------------------------------------------------------------------
 operation +proj=merc   +ellps=GRS80
 -------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 0 m
+accept  0 0
+expect  0 0
+tolerance 50 nm
 accept  2 1
-expect  222638.981586547 110579.965218250
+expect  222638.981586547 110579.965218249
 accept  2 -1
 expect  222638.981586547 -110579.965218249
 accept  -2 1
-expect  -222638.981586547 110579.965218250
+expect  -222638.981586547 110579.965218249
 accept  -2 -1
 expect  -222638.981586547 -110579.965218249
+# inflate tolerance by scale (k = 5.7e15)
+tolerance 3e8
+accept  0 89.99999999999999
+expect  0 235805185.015130176
+accept  0 -89.99999999999999
+expect  0 -235805185.015130176
 
 direction inverse
+tolerance 0 m
+accept  0 0
+expect  0 0
+tolerance 50 nm
 accept  200 100
-expect  0.001796631 0.000904369
+expect  0.00179663056824 0.00090436947704
 accept  200 -100
-expect  0.001796631 -0.000904369
+expect  0.00179663056824 -0.00090436947704
 accept  -200 100
-expect  -0.001796631 0.000904369
+expect  -0.00179663056824 0.00090436947704
 accept  -200 -100
-expect  -0.001796631 -0.000904369
+expect  -0.00179663056824 -0.00090436947704
+accept  0 235805185.015130176
+expect  0 89.99999999999999
+accept  0 -235805185.015130176
+expect  0 -89.99999999999999
+accept  0 1e10
+expect  0 90
+accept  0 -1e10
+expect  0 -90
 
 -------------------------------------------------------------------------------
 operation +proj=merc   +R=6400000
 -------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 0 m
+accept  0 0
+expect  0 0
+tolerance 50 nm
 accept  2 1
 expect  223402.144255274 111706.743574944
 accept  2 -1
@@ -3389,25 +3424,32 @@ accept  -2 -1
 expect  -223402.144255274 -111706.743574944
 
 direction inverse
+tolerance 0 m
+accept  0 0
+expect  0 0
+tolerance 50 nm
 accept  200 100
-expect  0.001790493 0.000895247
+expect  0.00179049310978 0.00089524655486
 accept  200 -100
-expect  0.001790493 -0.000895247
+expect  0.00179049310978 -0.00089524655486
 accept  -200 100
-expect  -0.001790493 0.000895247
+expect  -0.00179049310978 0.00089524655486
 accept  -200 -100
-expect  -0.001790493 -0.000895247
-
+expect  -0.00179049310978 -0.00089524655486
 
 -------------------------------------------------------------------------------
 operation +proj=merc +R=1
 -------------------------------------------------------------------------------
 # Test the numerical stability of the inverse spherical Mercator
 -------------------------------------------------------------------------------
-tolerance   1e-15 m
-accept  0   1e-15
+tolerance   1e-17 m
+accept  0   57.295779513e-15
 expect  0   1e-15
 
+direction inverse
+accept  0   1e-15
+expect  0   57.295779513e-15
+
 
 ===============================================================================
 # Miller Oblated Stereographic
@@ -5377,6 +5419,32 @@ accept  -200 -100
 expect  -0.001790493 -0.000895247
 
 
+-------------------------------------------------------------------------------
+operation +proj=stere   +ellps=GRS80        +lat_0=90  +lat_ts=70
+-------------------------------------------------------------------------------
+tolerance 1e-15m
+accept    0 90
+expect    0 0
+roundtrip 1
+
+tolerance 0.1 mm
+accept    20            70
+expect    748315.3282   -2055979.4669
+roundtrip 1
+
+-------------------------------------------------------------------------------
+operation +proj=stere   +ellps=GRS80        +lat_0=-90  +lat_ts=-70
+-------------------------------------------------------------------------------
+tolerance 1e-15m
+accept    0 -90
+expect    0 0
+roundtrip 1
+
+tolerance 0.1 mm
+accept    20            -70
+expect    748315.3282   2055979.4669
+roundtrip 1
+
 ===============================================================================
 # Oblique Stereographic Alternative
 # 	Azimuthal, Sph&Ell
@@ -5658,25 +5726,37 @@ expect  -0.001790143 0.511651393
 -------------------------------------------------------------------------------
 operation +proj=tmerc   +ellps=GRS80
 -------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 50 nm
 accept  2 1
-expect  222650.796795778 110642.229411927
+expect  222650.796797586 110642.229411933
 accept  2 -1
-expect  222650.796795778 -110642.229411927
+expect  222650.796797586 -110642.229411933
 accept  -2 1
-expect  -222650.796795778 110642.229411927
+expect  -222650.796797586 110642.229411933
 accept  -2 -1
-expect  -222650.796795778 -110642.229411927
+expect  -222650.796797586 -110642.229411933
+# near pole
+accept  30 89.9999
+expect  5.584698978 10001956.056248082
+# 3900 km from central meridian
+accept  44.69 35.37
+expect  4168136.489446198 4985511.302287407
 
 direction inverse
 accept  200 100
-expect  0.001796631 0.000904369
+expect  0.00179663056816 0.00090436947663
 accept  200 -100
-expect  0.001796631 -0.000904369
+expect  0.00179663056816 -0.00090436947663
 accept  -200 100
-expect  -0.001796631 0.000904369
+expect  -0.00179663056816 0.00090436947663
 accept  -200 -100
-expect  -0.001796631 -0.000904369
+expect  -0.00179663056816 -0.00090436947663
+# near pole
+accept  6 1.0001e7
+expect  0.35596960759234 89.99135362646302
+# 3900 km from central meridian
+accept  4168136.489446198 4985511.302287407
+expect  44.69 35.37
 
 -------------------------------------------------------------------------------
 operation +proj=tmerc   +R=6400000
@@ -5691,6 +5771,74 @@ expect  -223413.466406322 111769.145040597
 accept  -2 -1
 expect  -223413.466406322 -111769.145040597
 
+accept    0      0
+expect    0      0
+roundtrip 1
+
+accept    0      90
+expect    0      10053096.491487337276
+roundtrip 1
+
+accept    0      -90
+expect    0      -10053096.491487337276
+roundtrip 1
+
+accept    170                       60
+expect    557076.820490954560       13361866.764138307422
+roundtrip 1
+
+accept    89                        0.01
+expect    30344312.098578717560     64001.116414904580
+roundtrip 1
+
+accept    91                        0.01
+expect    30344312.098578717560     20042191.866555366665
+roundtrip 1
+
+accept    179.999                   0.01
+expect    111.701070433669          20105075.972255337983
+roundtrip 1
+
+accept    -89                       0.01
+expect    -30344312.098578717560    64001.116414904580
+roundtrip 1
+
+accept    -91                       0.01
+expect    -30344312.098578717560    20042191.866555366665
+roundtrip 1
+
+accept    -179.999                  0.01
+expect    -111.701070433669         20105075.972255337983
+roundtrip 1
+
+accept    89                       -0.01
+expect    30344312.098578717560    -64001.116414904580
+roundtrip 1
+
+accept    91                       -0.01
+expect    30344312.098578717560    -20042191.866555366665
+roundtrip 1
+
+accept    179.999                  -0.01
+expect    111.701070433669         -20105075.972255337983
+roundtrip 1
+
+accept    -89                      -0.01
+expect    -30344312.098578717560   -64001.116414904580
+roundtrip 1
+
+accept    -91                       -0.01
+expect    -30344312.098578717560    -20042191.866555366665
+roundtrip 1
+
+accept    -179.999                  -0.01
+expect    -111.701070433669         -20105075.972255337983
+roundtrip 1
+
+accept    150                       0
+expect    3515559.323737951461      20106192.982974674553
+roundtrip 1
+
 direction inverse
 accept  200 100
 expect  0.001790493 0.000895247
@@ -5702,6 +5850,130 @@ accept  -200 -100
 expect  -0.001790493 -0.000895247
 
 -------------------------------------------------------------------------------
+# Approx tmerc on a almost spherical ellipsoid, lat_0 north hemisphere
+
+operation +proj=tmerc   +a=6400000 +rf=1e12 +k=0.9 +lat_0=40 +approx
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept    0             -30
+expect    0             -7037167.5440
+roundtrip 1
+
+accept    1             -30
+expect    87064.5795    -7037547.4590
+roundtrip 1
+
+accept    -1            -30
+expect    -87064.5795   -7037547.4590
+roundtrip 1
+
+accept    0             30
+expect    0             -1005309.6491
+roundtrip 1
+
+accept    0             40
+expect    0             0
+roundtrip 1
+
+accept    1             41
+expect    75872.2182    100965.3718
+roundtrip 1
+
+-------------------------------------------------------------------------------
+# Approx tmerc on a sphere, lat_0 north hemisphere
+
+operation +proj=tmerc   +R=6400000 +k=0.9 +lat_0=40
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept    0             -30
+expect    0             -7037167.5440
+roundtrip 1
+
+accept    1             -30
+expect    87064.5795    -7037547.4590
+roundtrip 1
+
+accept    -1            -30
+expect    -87064.5795   -7037547.4590
+roundtrip 1
+
+accept    0             30
+expect    0             -1005309.6491
+roundtrip 1
+
+accept    0             40
+expect    0             0
+roundtrip 1
+
+accept    1             41
+expect    75872.2182    100965.3718
+roundtrip 1
+
+-------------------------------------------------------------------------------
+# Approx tmerc on a almost spherical ellipsoid, lat_0 south hemisphere
+
+operation +proj=tmerc   +a=6400000 +rf=1e12 +k=0.9 +lat_0=-40 +approx
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept    0             -30
+expect    0             1005309.6491
+roundtrip 1
+
+accept    1             -30
+expect    87064.5795    1004929.7341
+roundtrip 1
+
+accept    -1            -30
+expect    -87064.5795   1004929.7341
+roundtrip 1
+
+accept    0             30
+expect    0             7037167.5440
+roundtrip 1
+
+accept    0             -40
+expect    0             0
+roundtrip 1
+
+accept    1             -41
+expect    75872.2182    -100965.3718
+roundtrip 1
+
+-------------------------------------------------------------------------------
+# Approx tmerc on a sphere, lat_0 south hemisphere
+
+operation +proj=tmerc   +R=6400000 +k=0.9 +lat_0=-40
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept    0             -30
+expect    0             1005309.6491
+roundtrip 1
+
+accept    1             -30
+expect    87064.5795    1004929.7341
+roundtrip 1
+
+accept    -1            -30
+expect    -87064.5795   1004929.7341
+roundtrip 1
+
+accept    0             30
+expect    0             7037167.5440
+roundtrip 1
+
+accept    0             -40
+expect    0             0
+roundtrip 1
+
+accept    1             -41
+expect    75872.2182    -100965.3718
+roundtrip 1
+
+-------------------------------------------------------------------------------
 operation +proj=tmerc   +R=1
 -------------------------------------------------------------------------------
 direction inverse
@@ -6601,4 +6873,46 @@ accept    -74.25    4.8
 expect    80859.033 122543.174
 roundtrip 1
 
+===============================================================================
+# Geocentric/topocentric conversion
+===============================================================================
+
+# Test parameter and point from IOGP Publication 373-7-2 - Geomatics Guidance Note number 7, part 2 - October 2020
+operation +proj=topocentric +ellps=WGS84 +X_0=3652755.3058 +Y_0=319574.6799 +Z_0=5201547.3536
+tolerance 1 mm
+accept    3771793.968  140253.342         5124304.349
+expect    -189013.869  -128642.040        -4220.171
+roundtrip 1
+
+===============================================================================
+# Geographic/topocentric conversion
+===============================================================================
+
+# Test parameter and point from IOGP Publication 373-7-2 - Geomatics Guidance Note number 7, part 2 - October 2020
+operation +proj=pipeline +step +proj=cart +ellps=WGS84 +step +proj=topocentric +ellps=WGS84 +lon_0=5 +lat_0=55 +h_0=200
+tolerance 1 mm
+accept    2.12955      53.80939444444444  73
+expect    -189013.869  -128642.040        -4220.171
+roundtrip 1
+
+===============================================================================
+# Error cases of topocentric
+===============================================================================
+
+# missing X_0,Y_0,Z_0 or lon_0,lat_0
+operation +proj=topocentric +ellps=WGS84
+expect failure errno missing_args
+
+# missing Z_0
+operation +proj=topocentric +ellps=WGS84 +X_0=0 +Y_0=0
+expect failure errno missing_args
+
+# missing lat_0
+operation +proj=topocentric +ellps=WGS84 +lon_0=0
+expect failure errno missing_args
+
+# X_0 and lon_0 are mutually exclusive
+operation +proj=topocentric +ellps=WGS84 +X_0=0 +lon_0=0
+expect failure errno mutually_exclusive_args
+
 </gie-strict>
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 2c0c19a9..3924f47d 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -69,17 +69,6 @@ set(PROJ_TEST_ENVIRONMENT
   "PROJ_SOURCE_DATA=${PROJ_SOURCE_DIR}/data"
 )
 
-add_executable(proj_pj_transform_test
-  main.cpp
-  pj_transform_test.cpp)
-target_link_libraries(proj_pj_transform_test
-  GTest::gtest
-  ${PROJ_LIBRARIES})
-add_test(NAME proj_pj_transform_test COMMAND proj_pj_transform_test)
-set_property(TEST proj_pj_transform_test
-  PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT})
-
-
 add_executable(proj_errno_string_test
   main.cpp
   proj_errno_string_test.cpp)
diff --git a/test/unit/Makefile.am b/test/unit/Makefile.am
index 483cb0bd..4e931c2b 100644
--- a/test/unit/Makefile.am
+++ b/test/unit/Makefile.am
@@ -9,8 +9,7 @@ AM_CXXFLAGS = @CXX_WFLAGS@ @NO_ZERO_AS_NULL_POINTER_CONSTANT_FLAG@
 
 PROJ_LIB ?= ../../data/for_tests
 
-noinst_PROGRAMS = pj_transform_test
-noinst_PROGRAMS += pj_phi2_test
+noinst_PROGRAMS = pj_phi2_test
 noinst_PROGRAMS += proj_errno_string_test
 noinst_PROGRAMS += proj_angular_io_test
 noinst_PROGRAMS += proj_context_test
@@ -21,11 +20,6 @@ noinst_PROGRAMS += test_network
 noinst_PROGRAMS += test_defmodel
 noinst_PROGRAMS += test_tinshift
 
-pj_transform_test_SOURCES = pj_transform_test.cpp main.cpp
-pj_transform_test_LDADD = ../../src/libproj.la @GTEST_LIBS@
-
-pj_transform_test-check: pj_transform_test
-	PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) ./pj_transform_test
 
 pj_phi2_test_SOURCES = pj_phi2_test.cpp main.cpp
 pj_phi2_test_LDADD = ../../src/libproj.la @GTEST_LIBS@
@@ -84,6 +78,6 @@ test_tinshift_LDADD = ../../src/libproj.la @GTEST_LIBS@
 test_tinshift-check: test_tinshift
 	PROJ_LIB=$(PROJ_LIB) ./test_tinshift
 
-check-local: pj_transform_test-check pj_phi2_test-check proj_errno_string_test-check \
+check-local: pj_phi2_test-check proj_errno_string_test-check \
 		proj_angular_io_test-check proj_context_test-check test_cpp_api-check \
 		gie_self_tests-check test_network-check test_defmodel-check test_tinshift-check
diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
index 6f1b3c32..f9252137 100644
--- a/test/unit/gie_self_tests.cpp
+++ b/test/unit/gie_self_tests.cpp
@@ -548,7 +548,7 @@ static void test_time(const char *args, double tol, double t_in, double t_exp) {
     out = proj_trans(P, PJ_INV, out);
     EXPECT_NEAR(out.xyzt.t, t_in, tol);
 
-    pj_free(P);
+    proj_destroy(P);
 
     proj_log_level(NULL, PJ_LOG_NONE);
 }
diff --git a/test/unit/pj_phi2_test.cpp b/test/unit/pj_phi2_test.cpp
index c4db6e52..7ccbb01c 100644
--- a/test/unit/pj_phi2_test.cpp
+++ b/test/unit/pj_phi2_test.cpp
@@ -37,49 +37,58 @@
 namespace {
 
 TEST(PjPhi2Test, Basic) {
-    projCtx ctx = pj_get_default_ctx();
+    PJ_CONTEXT *ctx = pj_get_default_ctx();
+
+    // Expectation is that only sane values of e (and nan is here reckoned to
+    // be sane) are passed to pj_phi2.  Thus the return value with other values
+    // of e is "implementation dependent".
+
+    constexpr auto inf = std::numeric_limits<double>::infinity();
+    constexpr auto nan = std::numeric_limits<double>::quiet_NaN();
+
+    // Strict equality is demanded here.
+    EXPECT_EQ(M_PI_2, pj_phi2(ctx, +0.0, 0.0));
+    EXPECT_EQ(0.0, pj_phi2(ctx, 1.0, 0.0));
+    EXPECT_EQ(-M_PI_2, pj_phi2(ctx, inf, 0.0));
+    // We don't expect pj_phi2 to be called with negative ts (since ts =
+    // exp(-psi)).  However, in the current implementation it is odd in ts.
+    // N.B. ts = +0.0 and ts = -0.0 return different results.
+    EXPECT_EQ(-M_PI_2, pj_phi2(ctx, -0.0, 0.0));
+    EXPECT_EQ(0.0, pj_phi2(ctx, -1.0, 0.0));
+    EXPECT_EQ(+M_PI_2, pj_phi2(ctx, -inf, 0.0));
+
+    constexpr double e = 0.2;
+    EXPECT_EQ(M_PI_2, pj_phi2(ctx, +0.0, e));
+    EXPECT_EQ(0.0, pj_phi2(ctx, 1.0, e));
+    EXPECT_EQ(-M_PI_2, pj_phi2(ctx, inf, e));
+    EXPECT_EQ(-M_PI_2, pj_phi2(ctx, -0.0, e));
+    EXPECT_EQ(0.0, pj_phi2(ctx, -1.0, e));
+    EXPECT_EQ(+M_PI_2, pj_phi2(ctx, -inf, e));
 
-    EXPECT_DOUBLE_EQ(M_PI_2, pj_phi2(ctx, 0.0, 0.0));
-
-    EXPECT_NEAR(0.0, pj_phi2(ctx, 1.0, 0.0), 1e-16);
-    EXPECT_DOUBLE_EQ(M_PI_2, pj_phi2(ctx, 0.0, 1.0));
-    EXPECT_DOUBLE_EQ(M_PI, pj_phi2(ctx, -1.0, 0.0));
-    EXPECT_DOUBLE_EQ(M_PI_2, pj_phi2(ctx, 0.0, -1.0));
-
-    EXPECT_NEAR(0.0, pj_phi2(ctx, 1.0, 1.0), 1e-16);
-    EXPECT_DOUBLE_EQ(M_PI, pj_phi2(ctx, -1.0, -1.0));
-
-    // TODO(schwehr): M_PI_4, M_PI_2, M_PI, M_E
-    // https://www.gnu.org/software/libc/manual/html_node/Mathematical-Constants.html
-
-    EXPECT_DOUBLE_EQ(-0.95445818456292697, pj_phi2(ctx, M_PI, 0.0));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, M_PI)));
-    EXPECT_DOUBLE_EQ(4.0960508381527205, pj_phi2(ctx, -M_PI, 0.0));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, -M_PI)));
-
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, M_PI, M_PI)));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -M_PI, -M_PI)));
-}
-
-TEST(PjPhi2Test, AvoidUndefinedBehavior) {
-    auto ctx = pj_get_default_ctx();
-
-    const auto nan = std::numeric_limits<double>::quiet_NaN();
     EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, 0.0)));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, e)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, +0.0, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 1.0, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, inf, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -0.0, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -1.0, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -inf, nan)));
     EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, nan)));
 
-    // We do not really care about the values that follow.
-    const auto inf = std::numeric_limits<double>::infinity();
-
-    EXPECT_DOUBLE_EQ(-M_PI_2, pj_phi2(ctx, inf, 0.0));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, inf)));
-
-    EXPECT_DOUBLE_EQ(4.7123889803846897, pj_phi2(ctx, -inf, 0.0));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, -inf)));
-
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, inf, inf)));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -inf, -inf)));
+    EXPECT_DOUBLE_EQ(M_PI / 3, pj_phi2(ctx, 1 / (sqrt(3.0) + 2), 0.0));
+    EXPECT_DOUBLE_EQ(M_PI / 4, pj_phi2(ctx, 1 / (sqrt(2.0) + 1), 0.0));
+    EXPECT_DOUBLE_EQ(M_PI / 6, pj_phi2(ctx, 1 / sqrt(3.0), 0.0));
+    EXPECT_DOUBLE_EQ(-M_PI / 3, pj_phi2(ctx, sqrt(3.0) + 2, 0.0));
+    EXPECT_DOUBLE_EQ(-M_PI / 4, pj_phi2(ctx, sqrt(2.0) + 1, 0.0));
+    EXPECT_DOUBLE_EQ(-M_PI / 6, pj_phi2(ctx, sqrt(3.0), 0.0));
+
+    // Generated with exp(e * atanh(e * sin(phi))) / (tan(phi) + sec(phi))
+    EXPECT_DOUBLE_EQ(M_PI / 3, pj_phi2(ctx, 0.27749174377027023413, e));
+    EXPECT_DOUBLE_EQ(M_PI / 4, pj_phi2(ctx, 0.42617788119104192995, e));
+    EXPECT_DOUBLE_EQ(M_PI / 6, pj_phi2(ctx, 0.58905302448626726064, e));
+    EXPECT_DOUBLE_EQ(-M_PI / 3, pj_phi2(ctx, 3.6037108218537833089, e));
+    EXPECT_DOUBLE_EQ(-M_PI / 4, pj_phi2(ctx, 2.3464380582241712935, e));
+    EXPECT_DOUBLE_EQ(-M_PI / 6, pj_phi2(ctx, 1.6976400399134411849, e));
 }
 
 } // namespace
diff --git a/test/unit/pj_transform_test.cpp b/test/unit/pj_transform_test.cpp
deleted file mode 100644
index ddb054f0..00000000
--- a/test/unit/pj_transform_test.cpp
+++ /dev/null
@@ -1,740 +0,0 @@
-/******************************************************************************
- *
- * Project:  PROJ
- * Purpose:  Test pj_transform() legacy interface
- * Author:   Even Rouault <even dot rouault at spatialys dot com>
- *
- ******************************************************************************
- * Copyright (c) 2018, Even Rouault <even dot rouault at spatialys dot com>
- *
- * 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.
- ****************************************************************************/
-
-#define ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
-
-#include "gtest_include.h"
-#include <memory>
-
-// PROJ include order is sensitive
-// clang-format off
-#include <proj.h>
-#include "proj_internal.h"
-#include <proj_api.h>
-// clang-format on
-
-namespace {
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_to_longlat) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    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_EQ(x, 2 * DEG_TO_RAD);
-    EXPECT_EQ(y, 49 * DEG_TO_RAD);
-
-    x = 182 * DEG_TO_RAD;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
-    EXPECT_EQ(x, 182 * DEG_TO_RAD);
-    EXPECT_EQ(y, 49 * DEG_TO_RAD);
-
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_to_proj) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst = pj_init_plus("+proj=utm +zone=31 +datum=WGS84");
-    double x = 3 * DEG_TO_RAD;
-    double y = 0 * DEG_TO_RAD;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
-    EXPECT_NEAR(x, 500000, 1e-8);
-    EXPECT_NEAR(y, 0, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_to_proj_tometer) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst = pj_init_plus("+proj=utm +zone=31 +datum=WGS84 +to_meter=1000");
-    double x = 3 * DEG_TO_RAD;
-    double y = 0 * DEG_TO_RAD;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
-    EXPECT_NEAR(x, 500, 1e-8);
-    EXPECT_NEAR(y, 0, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, proj_to_longlat) {
-    auto src = pj_init_plus("+proj=utm +zone=31 +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    double x = 500000;
-    double y = 0;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
-    EXPECT_NEAR(x, 3 * DEG_TO_RAD, 1e-12);
-    EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, proj_to_proj) {
-    auto src = pj_init_plus("+proj=utm +zone=31 +datum=WGS84");
-    auto dst = pj_init_plus("+proj=utm +zone=31 +datum=WGS84");
-    double x = 500000;
-    double y = 0;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
-    EXPECT_NEAR(x, 500000, 1e-8);
-    EXPECT_NEAR(y, 0, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_to_geocent) {
-    auto src = pj_init_plus("+proj=longlat +R=2");
-    auto dst = pj_init_plus("+proj=geocent +R=2");
-    double x = 0;
-    double y = 0;
-    double z = 0;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 2, 1e-8);
-    EXPECT_NEAR(y, 0, 1e-8);
-    EXPECT_NEAR(z, 0, 1e-8);
-
-    x = 90 * DEG_TO_RAD;
-    y = 0;
-    z = 0;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 0, 1e-8);
-    EXPECT_NEAR(y, 2, 1e-8);
-    EXPECT_NEAR(z, 0, 1e-8);
-
-    x = 0;
-    y = 90 * DEG_TO_RAD;
-    z = 0.1;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 0, 1e-8);
-    EXPECT_NEAR(y, 0, 1e-8);
-    EXPECT_NEAR(z, 2 + 0.1, 1e-8);
-
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_to_geocent_to_meter) {
-    auto src = pj_init_plus("+proj=longlat +R=2");
-    auto dst = pj_init_plus("+proj=geocent +R=2 +to_meter=1000");
-    double x = 0;
-    double y = 0;
-    double z = 0;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 2e-3, 1e-8);
-    EXPECT_NEAR(y, 0, 1e-8);
-    EXPECT_NEAR(z, 0, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, geocent_to_longlat) {
-    auto src = pj_init_plus("+proj=geocent +R=2");
-    auto dst = pj_init_plus("+proj=longlat +R=2");
-    double x = 0;
-    double y = 2;
-    double z = 0;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 0, 1e-12) << y / DEG_TO_RAD;
-    EXPECT_NEAR(z, 0, 1e-12);
-
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, geocent_to_meter_to_longlat) {
-    auto src = pj_init_plus("+proj=geocent +to_meter=1000 +R=2");
-    auto dst = pj_init_plus("+proj=longlat +R=2");
-    double x = 0;
-    double y = 2e-3;
-    double z = 0;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 0, 1e-12) << y / DEG_TO_RAD;
-    EXPECT_NEAR(z, 0, 1e-12);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, pm) {
-    auto src = pj_init_plus("+proj=longlat +pm=3 +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +pm=1 +datum=WGS84");
-    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 + 3 - 1) * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_EQ(y, 49 * DEG_TO_RAD) << y / DEG_TO_RAD;
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_geoc_to_longlat) {
-    auto src = pj_init_plus("+proj=longlat +geoc +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    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, 48.809360314691766 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_to_longlat_geoc) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +geoc +datum=WGS84");
-    double x = 2 * DEG_TO_RAD;
-    double y = 48.809360314691766 * 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);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, ellps_to_ellps_noop) {
-    auto src = pj_init_plus("+proj=longlat +ellps=clrk66");
-    auto dst = pj_init_plus("+proj=longlat +ellps=WGS84");
-    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);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, towgs84_3param_noop) {
-    auto src = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=1,2,3");
-    auto dst = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=1,2,3");
-    double x = 2 * DEG_TO_RAD;
-    double y = 49 * DEG_TO_RAD;
-    double z = 10;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 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;
-    EXPECT_NEAR(z, 10, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, towgs84_7param_noop) {
-    auto src =
-        pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=1,2,3,4,5,6,7");
-    auto dst =
-        pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=1,2,3,4,5,6,7");
-    double x = 2 * DEG_TO_RAD;
-    double y = 49 * DEG_TO_RAD;
-    double z = 10;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 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;
-    EXPECT_NEAR(z, 10, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_towgs84_3param_to_datum) {
-    auto src = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    double x = 90 * DEG_TO_RAD;
-    double y = 0 * DEG_TO_RAD;
-    double z = 10;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-    EXPECT_NEAR(z, 10 + 1, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_towgs84_3param_to_datum_no_z) {
-    auto src = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    double x = 90 * DEG_TO_RAD;
-    double y = 0 * DEG_TO_RAD;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
-    EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_towgs84_7param_to_datum) {
-    auto src =
-        pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0,0,0,0,0.5");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    double x = 90 * DEG_TO_RAD;
-    double y = 0 * DEG_TO_RAD;
-    double z = 10;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-    EXPECT_NEAR(z, 14.189073500223458, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, datum_to_longlat_towgs84_3param) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0");
-    double x = 90 * DEG_TO_RAD;
-    double y = 0 * DEG_TO_RAD;
-    double z = 10 + 1;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-    EXPECT_NEAR(z, 10, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, datum_to_longlat_towgs84_7param) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst =
-        pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0,0,0,0,0.5");
-    double x = 90 * DEG_TO_RAD;
-    double y = 0 * DEG_TO_RAD;
-    double z = 14.189073500223458;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-    EXPECT_NEAR(z, 10, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, ellps_grs80_towgs84_to_datum_wgs84) {
-    auto src = pj_init_plus("+proj=longlat +ellps=GRS80 +towgs84=0,0,0");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    double x = 2 * DEG_TO_RAD;
-    double y = 49 * DEG_TO_RAD;
-    double z = 10;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-15) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-15) << y / DEG_TO_RAD;
-    EXPECT_NEAR(z, 10, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, longlat_nadgrids_to_datum) {
-    auto src = pj_init_plus("+proj=longlat +ellps=clrk66 +nadgrids=conus");
-    auto dst = pj_init_plus("+proj=longlat +datum=NAD83");
-    double x = -100 * DEG_TO_RAD;
-    double y = 40 * DEG_TO_RAD;
-    double z = 10;
-    int ret = pj_transform(src, dst, 1, 0, &x, &y, &z);
-    EXPECT_TRUE(ret == 0 || ret == PJD_ERR_FAILED_TO_LOAD_GRID);
-    if (ret == 0) {
-        EXPECT_NEAR(x, -100.00040583667015 * DEG_TO_RAD, 1e-12)
-            << x / DEG_TO_RAD;
-        EXPECT_NEAR(y, 40.000005895651363 * DEG_TO_RAD, 1e-12)
-            << y / DEG_TO_RAD;
-        EXPECT_NEAR(z, 10.000043224543333, 1e-8);
-    }
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, nadgrids_noop) {
-    auto src = pj_init_plus("+proj=longlat +ellps=clrk66 +nadgrids=conus");
-    auto dst = pj_init_plus("+proj=longlat +ellps=clrk66 +nadgrids=conus");
-    double x = -100 * DEG_TO_RAD;
-    double y = 40 * DEG_TO_RAD;
-    double z = 10;
-    int ret = pj_transform(src, dst, 1, 0, &x, &y, &z);
-    EXPECT_TRUE(ret == 0);
-    if (ret == 0) {
-        EXPECT_NEAR(x, -100 * DEG_TO_RAD, 1e-15) << x / DEG_TO_RAD;
-        EXPECT_NEAR(y, 40 * DEG_TO_RAD, 1e-15) << y / DEG_TO_RAD;
-        EXPECT_NEAR(z, 10, 1e-8);
-    }
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, datum_to_longlat_nadgrids) {
-    auto src = pj_init_plus("+proj=longlat +datum=NAD83");
-    auto dst = pj_init_plus("+proj=longlat +ellps=clrk66 +nadgrids=conus");
-    double x = -100.00040583667015 * DEG_TO_RAD;
-    double y = 40.000005895651363 * DEG_TO_RAD;
-    double z = 10.000043224543333;
-    int ret = pj_transform(src, dst, 1, 0, &x, &y, &z);
-    EXPECT_TRUE(ret == 0 || ret == PJD_ERR_FAILED_TO_LOAD_GRID);
-    if (ret == 0) {
-        EXPECT_NEAR(x, -100 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-        EXPECT_NEAR(y, 40 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-        EXPECT_NEAR(z, 10, 1e-8);
-    }
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, long_wrap) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84 +lon_wrap=180");
-    double x = -1 * DEG_TO_RAD;
-    double y = 0 * DEG_TO_RAD;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
-    EXPECT_NEAR(x, 359 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, src_vto_meter) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84 +vto_meter=1000");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    double x = 2 * DEG_TO_RAD;
-    double y = 49 * DEG_TO_RAD;
-    double z = 1;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 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;
-    EXPECT_NEAR(z, 1000, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, dest_vto_meter) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84 +vto_meter=1000");
-    double x = 2 * DEG_TO_RAD;
-    double y = 49 * DEG_TO_RAD;
-    double z = 1000;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 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;
-    EXPECT_NEAR(z, 1, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, src_axis_neu_with_z) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84 +axis=neu");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    double x = 49 * DEG_TO_RAD;
-    double y = 2 * DEG_TO_RAD;
-    double z = 1;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 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;
-    EXPECT_NEAR(z, 1, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, src_axis_neu_without_z) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84 +axis=neu");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    double x = 49 * DEG_TO_RAD;
-    double y = 2 * 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);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, src_axis_swd) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84 +axis=swd");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84");
-    double x = 49 * DEG_TO_RAD;
-    double y = 2 * DEG_TO_RAD;
-    double z = -1;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 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;
-    EXPECT_NEAR(z, 1, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, dst_axis_neu) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84 +axis=neu");
-    double x = 2 * DEG_TO_RAD;
-    double y = 49 * DEG_TO_RAD;
-    double z = 1;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, 49 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, 2 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-    EXPECT_NEAR(z, 1, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, dst_axis_swd) {
-    auto src = pj_init_plus("+proj=longlat +datum=WGS84");
-    auto dst = pj_init_plus("+proj=longlat +datum=WGS84 +axis=swd");
-    double x = 2 * DEG_TO_RAD;
-    double y = 49 * DEG_TO_RAD;
-    double z = 1;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0);
-    EXPECT_NEAR(x, -49 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD;
-    EXPECT_NEAR(y, -2 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD;
-    EXPECT_NEAR(z, -1, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(pj_transform_test, init_epsg) {
-    auto src = pj_init_plus("+init=epsg:4326");
-    ASSERT_TRUE(src != nullptr);
-    auto dst = pj_init_plus("+init=epsg:32631");
-    ASSERT_TRUE(dst != nullptr);
-    double x = 3 * DEG_TO_RAD;
-    double y = 0 * DEG_TO_RAD;
-    EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0);
-    EXPECT_NEAR(x, 500000, 1e-8);
-    EXPECT_NEAR(y, 0, 1e-8);
-    pj_free(src);
-    pj_free(dst);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(proj_api_h, pj_set_searchpath) {
-
-    const char *path = "/i_do/not/exit";
-    pj_set_searchpath(1, &path);
-    {
-        auto info = proj_info();
-        EXPECT_EQ(info.path_count, 1U);
-        ASSERT_NE(info.paths, nullptr);
-        ASSERT_NE(info.paths[0], nullptr);
-        EXPECT_EQ(std::string(info.paths[0]), path);
-    }
-
-    pj_set_searchpath(0, nullptr);
-    {
-        auto info = proj_info();
-        EXPECT_EQ(info.path_count, 0U);
-        EXPECT_EQ(info.paths, nullptr);
-    }
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(proj_api_h, pj_set_finder) {
-
-    const auto myfinder = [](const char *) -> const char * { return nullptr; };
-    pj_set_finder(myfinder);
-
-    pj_set_finder(nullptr);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(proj_api_h, default_fileapi) {
-    auto ctx = pj_ctx_alloc();
-    auto fid = pj_open_lib(ctx, "proj.db", "rb");
-    ASSERT_NE(fid, nullptr);
-    char header[6];
-    ASSERT_EQ(pj_ctx_fread(ctx, header, 1, 6, fid), 6U);
-    ASSERT_TRUE(memcmp(header, "SQLite", 6) == 0);
-    ASSERT_EQ(pj_ctx_ftell(ctx, fid), 6);
-    ASSERT_EQ(pj_ctx_fseek(ctx, fid, 0, SEEK_SET), 0);
-    ASSERT_EQ(pj_ctx_ftell(ctx, fid), 0);
-    pj_ctx_fclose(ctx, fid);
-    pj_ctx_free(ctx);
-}
-
-// ---------------------------------------------------------------------------
-
-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;
-    bool gotInMyFSeek = false;
-    bool gotInMyFTell = false;
-    bool gotInMyFClose = false;
-};
-
-struct MyFile {
-    FILE *fp;
-    Spy *spy;
-};
-
-static PAFile myFOpen(projCtx ctx, const char *filename, const char *access) {
-    FILE *fp = fopen(filename, access);
-    if (!fp)
-        return nullptr;
-    MyFile *myF = new MyFile;
-    myF->spy = (Spy *)pj_ctx_get_app_data(ctx);
-    myF->spy->gotInMyFOpen = true;
-    myF->fp = fp;
-    return reinterpret_cast<PAFile>(myF);
-}
-
-static size_t myFRead(void *buffer, size_t size, size_t nmemb, PAFile file) {
-    MyFile *myF = reinterpret_cast<MyFile *>(file);
-    myF->spy->gotInMyFRead = true;
-    return fread(buffer, size, nmemb, myF->fp);
-}
-
-static int myFSeek(PAFile file, long offset, int whence) {
-    MyFile *myF = reinterpret_cast<MyFile *>(file);
-    myF->spy->gotInMyFSeek = true;
-    return fseek(myF->fp, offset, whence);
-}
-
-static long myFTell(PAFile file) {
-    MyFile *myF = reinterpret_cast<MyFile *>(file);
-    myF->spy->gotInMyFTell = true;
-    return ftell(myF->fp);
-}
-
-static void myFClose(PAFile file) {
-    MyFile *myF = reinterpret_cast<MyFile *>(file);
-    myF->spy->gotInMyFClose = true;
-    fclose(myF->fp);
-    delete myF;
-}
-
-TEST(proj_api_h, custom_fileapi) {
-    auto ctx = pj_ctx_alloc();
-    Spy spy;
-    pj_ctx_set_app_data(ctx, &spy);
-    projFileAPI myAPI = {myFOpen, myFRead, myFSeek, myFTell, myFClose};
-    pj_ctx_set_fileapi(ctx, &myAPI);
-    EXPECT_EQ(pj_ctx_get_fileapi(ctx), &myAPI);
-    auto fid = pj_open_lib(ctx, "proj.db", "rb");
-    ASSERT_NE(fid, nullptr);
-    char header[6];
-    ASSERT_EQ(pj_ctx_fread(ctx, header, 1, 6, fid), 6U);
-    ASSERT_TRUE(memcmp(header, "SQLite", 6) == 0);
-    ASSERT_EQ(pj_ctx_ftell(ctx, fid), 6);
-    ASSERT_EQ(pj_ctx_fseek(ctx, fid, 0, SEEK_SET), 0);
-    ASSERT_EQ(pj_ctx_ftell(ctx, fid), 0);
-    pj_ctx_fclose(ctx, fid);
-    pj_ctx_free(ctx);
-    EXPECT_TRUE(spy.gotInMyFOpen);
-    EXPECT_TRUE(spy.gotInMyFRead);
-    EXPECT_TRUE(spy.gotInMyFSeek);
-    EXPECT_TRUE(spy.gotInMyFTell);
-    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 5da6b369..c417371d 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -433,16 +433,27 @@ TEST_F(CApi, proj_as_wkt) {
     ObjectKeeper keeper_crs4979(crs4979);
     ASSERT_NE(crs4979, nullptr);
 
+    EXPECT_EQ(proj_as_wkt(m_ctxt, crs4979, PJ_WKT1_GDAL, nullptr), nullptr);
+
     // STRICT=NO
     {
-        EXPECT_EQ(proj_as_wkt(m_ctxt, crs4979, PJ_WKT1_GDAL, nullptr), nullptr);
-
         const char *const options[] = {"STRICT=NO", nullptr};
         auto wkt = proj_as_wkt(m_ctxt, crs4979, PJ_WKT1_GDAL, options);
         ASSERT_NE(wkt, nullptr);
         EXPECT_TRUE(std::string(wkt).find("GEOGCS[\"WGS 84\"") == 0) << wkt;
     }
 
+    // ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS=YES
+    {
+        const char *const options[] = {
+            "ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS=YES", nullptr};
+        auto wkt = proj_as_wkt(m_ctxt, crs4979, PJ_WKT1_GDAL, options);
+        ASSERT_NE(wkt, nullptr);
+        EXPECT_TRUE(std::string(wkt).find(
+                        "COMPD_CS[\"WGS 84 + Ellipsoid (metre)\"") == 0)
+            << wkt;
+    }
+
     // unsupported option
     {
         const char *const options[] = {"unsupported=yes", nullptr};
@@ -845,6 +856,13 @@ TEST_F(CApi, proj_create_from_database) {
         EXPECT_EQ(proj_get_type(datum), PJ_TYPE_GEODETIC_REFERENCE_FRAME);
     }
     {
+        auto ensemble = proj_create_from_database(
+            m_ctxt, "EPSG", "6326", PJ_CATEGORY_DATUM_ENSEMBLE, false, nullptr);
+        ASSERT_NE(ensemble, nullptr);
+        ObjectKeeper keeper(ensemble);
+        EXPECT_EQ(proj_get_type(ensemble), PJ_TYPE_DATUM_ENSEMBLE);
+    }
+    {
         // International Terrestrial Reference Frame 2008
         auto datum = proj_create_from_database(
             m_ctxt, "EPSG", "1061", PJ_CATEGORY_DATUM, false, nullptr);
@@ -4643,10 +4661,21 @@ TEST_F(CApi, proj_create_vertical_crs_ex) {
     ObjectKeeper keeper_geog_crs(geog_crs);
     ASSERT_NE(geog_crs, nullptr);
 
-    auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr,
-                                            nullptr);
-    ObjectKeeper keeper_P(P);
-    ASSERT_NE(P, nullptr);
+    PJ_OPERATION_FACTORY_CONTEXT *ctxt =
+        proj_create_operation_factory_context(m_ctxt, nullptr);
+    ASSERT_NE(ctxt, nullptr);
+    ContextKeeper keeper_ctxt(ctxt);
+    proj_operation_factory_context_set_grid_availability_use(
+        m_ctxt, ctxt, PROJ_GRID_AVAILABILITY_IGNORED);
+    proj_operation_factory_context_set_spatial_criterion(
+        m_ctxt, ctxt, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
+    PJ_OBJ_LIST *operations =
+        proj_create_operations(m_ctxt, compound, geog_crs, ctxt);
+    ASSERT_NE(operations, nullptr);
+    ObjListKeeper keeper_operations(operations);
+    EXPECT_GE(proj_list_get_count(operations), 1);
+    auto P = proj_list_get(m_ctxt, operations, 0);
+    ObjectKeeper keeper_transform(P);
 
     auto name = proj_get_name(P);
     ASSERT_TRUE(name != nullptr);
@@ -4699,10 +4728,21 @@ TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) {
     ObjectKeeper keeper_geog_crs(geog_crs);
     ASSERT_NE(geog_crs, nullptr);
 
-    auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr,
-                                            nullptr);
-    ObjectKeeper keeper_P(P);
-    ASSERT_NE(P, nullptr);
+    PJ_OPERATION_FACTORY_CONTEXT *ctxt =
+        proj_create_operation_factory_context(m_ctxt, nullptr);
+    ASSERT_NE(ctxt, nullptr);
+    ContextKeeper keeper_ctxt(ctxt);
+    proj_operation_factory_context_set_grid_availability_use(
+        m_ctxt, ctxt, PROJ_GRID_AVAILABILITY_IGNORED);
+    proj_operation_factory_context_set_spatial_criterion(
+        m_ctxt, ctxt, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
+    PJ_OBJ_LIST *operations =
+        proj_create_operations(m_ctxt, compound, geog_crs, ctxt);
+    ASSERT_NE(operations, nullptr);
+    ObjListKeeper keeper_operations(operations);
+    EXPECT_GE(proj_list_get_count(operations), 1);
+    auto P = proj_list_get(m_ctxt, operations, 0);
+    ObjectKeeper keeper_transform(P);
 
     auto name = proj_get_name(P);
     ASSERT_TRUE(name != nullptr);
@@ -4732,10 +4772,13 @@ TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) {
     ObjectKeeper keeper_compound_from_projjson(compound_from_projjson);
     ASSERT_NE(compound_from_projjson, nullptr);
 
-    auto P2 = proj_create_crs_to_crs_from_pj(m_ctxt, compound_from_projjson,
-                                             geog_crs, nullptr, nullptr);
-    ObjectKeeper keeper_P2(P2);
-    ASSERT_NE(P2, nullptr);
+    PJ_OBJ_LIST *operations2 =
+        proj_create_operations(m_ctxt, compound_from_projjson, geog_crs, ctxt);
+    ASSERT_NE(operations2, nullptr);
+    ObjListKeeper keeper_operations2(operations2);
+    EXPECT_GE(proj_list_get_count(operations2), 1);
+    auto P2 = proj_list_get(m_ctxt, operations2, 0);
+    ObjectKeeper keeper_transform2(P2);
 
     auto name_bis = proj_get_name(P2);
     ASSERT_TRUE(name_bis != nullptr);
@@ -4841,9 +4884,16 @@ TEST_F(CApi, proj_create_derived_geographic_crs) {
     const char *expected_wkt =
         "GEOGCRS[\"my rotated CRS\",\n"
         "    BASEGEOGCRS[\"WGS 84\",\n"
-        "        DATUM[\"World Geodetic System 1984\",\n"
+        "        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n"
+        "            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G730)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G873)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n"
         "            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
-        "                LENGTHUNIT[\"metre\",1]]],\n"
+        "                LENGTHUNIT[\"metre\",1]],\n"
+        "            ENSEMBLEACCURACY[2.0]],\n"
         "        PRIMEM[\"Greenwich\",0,\n"
         "            ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
         "    DERIVINGCONVERSION[\"Pole rotation (GRIB convention)\",\n"
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index a0fee905..cea5ff6c 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -414,6 +414,20 @@ TEST(crs, EPSG_4326_as_WKT1_ESRI_with_database) {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, EPSG_4901_as_WKT1_ESRI_with_PRIMEM_unit_name_morphing) {
+    auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto crs = factory->createCoordinateReferenceSystem("4901");
+    WKTFormatterNNPtr f(WKTFormatter::create(
+        WKTFormatter::Convention::WKT1_ESRI, DatabaseContext::create()));
+    EXPECT_EQ(crs->exportToWKT(f.get()),
+              "GEOGCS[\"GCS_ATF_Paris\",DATUM[\"D_ATF\","
+              "SPHEROID[\"Plessis_1817\",6376523.0,308.64]],"
+              "PRIMEM[\"Paris_RGS\",2.33720833333333],"
+              "UNIT[\"Grad\",0.0157079632679489]]");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, EPSG_4326_as_WKT1_ESRI_without_database) {
     auto crs = GeographicCRS::EPSG_4326;
     WKTFormatterNNPtr f(
@@ -513,6 +527,34 @@ TEST(crs, EPSG_4979_as_WKT1_GDAL) {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, EPSG_4979_as_WKT1_GDAL_with_ellipsoidal_height_as_vertical_crs) {
+    auto crs = GeographicCRS::EPSG_4979;
+    auto wkt = crs->exportToWKT(
+        &(WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL,
+                               DatabaseContext::create())
+              ->setAllowEllipsoidalHeightAsVerticalCRS(true)));
+
+    // For LAS 1.4 WKT1...
+    EXPECT_EQ(wkt, "COMPD_CS[\"WGS 84 + Ellipsoid (metre)\",\n"
+                   "    GEOGCS[\"WGS 84\",\n"
+                   "        DATUM[\"WGS_1984\",\n"
+                   "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+                   "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+                   "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+                   "        PRIMEM[\"Greenwich\",0,\n"
+                   "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+                   "        UNIT[\"degree\",0.0174532925199433,\n"
+                   "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+                   "        AUTHORITY[\"EPSG\",\"4326\"]],\n"
+                   "    VERT_CS[\"Ellipsoid (metre)\",\n"
+                   "        VERT_DATUM[\"Ellipsoid\",2002],\n"
+                   "        UNIT[\"metre\",1,\n"
+                   "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+                   "        AXIS[\"Ellipsoidal height\",UP]]]");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, EPSG_4979_as_WKT1_ESRI) {
     auto crs = GeographicCRS::EPSG_4979;
     WKTFormatterNNPtr f(
@@ -2037,6 +2079,50 @@ TEST(crs, projectedCRS_as_WKT1_ESRI) {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs,
+     projectedCRS_3D_as_WKT1_GDAL_with_ellipsoidal_height_as_vertical_crs) {
+    auto dbContext = DatabaseContext::create();
+    auto crs = AuthorityFactory::create(dbContext, "EPSG")
+                   ->createProjectedCRS("32631")
+                   ->promoteTo3D(std::string(), dbContext);
+    auto wkt = crs->exportToWKT(
+        &(WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext)
+              ->setAllowEllipsoidalHeightAsVerticalCRS(true)));
+
+    // For LAS 1.4 WKT1...
+    EXPECT_EQ(wkt,
+              "COMPD_CS[\"WGS 84 / UTM zone 31N + Ellipsoid (metre)\",\n"
+              "    PROJCS[\"WGS 84 / UTM zone 31N\",\n"
+              "        GEOGCS[\"WGS 84\",\n"
+              "            DATUM[\"WGS_1984\",\n"
+              "                SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+              "                    AUTHORITY[\"EPSG\",\"7030\"]],\n"
+              "                AUTHORITY[\"EPSG\",\"6326\"]],\n"
+              "            PRIMEM[\"Greenwich\",0,\n"
+              "                AUTHORITY[\"EPSG\",\"8901\"]],\n"
+              "            UNIT[\"degree\",0.0174532925199433,\n"
+              "                AUTHORITY[\"EPSG\",\"9122\"]],\n"
+              "            AUTHORITY[\"EPSG\",\"4326\"]],\n"
+              "        PROJECTION[\"Transverse_Mercator\"],\n"
+              "        PARAMETER[\"latitude_of_origin\",0],\n"
+              "        PARAMETER[\"central_meridian\",3],\n"
+              "        PARAMETER[\"scale_factor\",0.9996],\n"
+              "        PARAMETER[\"false_easting\",500000],\n"
+              "        PARAMETER[\"false_northing\",0],\n"
+              "        UNIT[\"metre\",1,\n"
+              "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+              "        AXIS[\"Easting\",EAST],\n"
+              "        AXIS[\"Northing\",NORTH],\n"
+              "        AUTHORITY[\"EPSG\",\"32631\"]],\n"
+              "    VERT_CS[\"Ellipsoid (metre)\",\n"
+              "        VERT_DATUM[\"Ellipsoid\",2002],\n"
+              "        UNIT[\"metre\",1,\n"
+              "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+              "        AXIS[\"Ellipsoidal height\",UP]]]");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, projectedCRS_with_ESRI_code_as_WKT1_ESRI) {
     auto dbContext = DatabaseContext::create();
     auto crs = AuthorityFactory::create(dbContext, "ESRI")
@@ -2107,6 +2193,53 @@ TEST(crs, projectedCRS_as_PROJ_string) {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, projectedCRS_3D_is_WKT1_equivalent_to_WKT2) {
+    auto dbContext = DatabaseContext::create();
+
+    // "Illegal" WKT1 with a Projected 3D CRS
+    auto wkt1 = "PROJCS[\"WGS 84 / UTM zone 16N [EGM08-1]\","
+                "GEOGCS[\"WGS 84 / UTM zone 16N [EGM08-1]\","
+                "DATUM[\"WGS84\",SPHEROID[\"WGS84\",6378137.000,298.257223563,"
+                "AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],"
+                "PRIMEM[\"Greenwich\",0.0000000000000000,"
+                "AUTHORITY[\"EPSG\",\"8901\"]],"
+                "UNIT[\"Degree\",0.01745329251994329547,"
+                "AUTHORITY[\"EPSG\",\"9102\"]],AUTHORITY[\"EPSG\",\"32616\"]],"
+                "UNIT[\"Meter\",1.00000000000000000000,"
+                "AUTHORITY[\"EPSG\",\"9001\"]],"
+                "PROJECTION[\"Transverse_Mercator\"],"
+                "PARAMETER[\"latitude_of_origin\",0.0000000000000000],"
+                "PARAMETER[\"central_meridian\",-87.0000000002777938],"
+                "PARAMETER[\"scale_factor\",0.9996000000000000],"
+                "PARAMETER[\"false_easting\",500000.000],"
+                "PARAMETER[\"false_northing\",0.000],"
+                "AXIS[\"Easting\",EAST],"
+                "AXIS[\"Northing\",NORTH],"
+                "AXIS[\"Height\",UP],"
+                "AUTHORITY[\"EPSG\",\"32616\"]]";
+
+    auto obj = WKTParser()
+                   .setStrict(false)
+                   .attachDatabaseContext(dbContext)
+                   .createFromWKT(wkt1);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    WKTFormatterNNPtr f(
+        WKTFormatter::create(WKTFormatter::Convention::WKT2_2019));
+    auto wkt2 = crs->exportToWKT(f.get());
+    auto obj2 =
+        WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt2);
+    auto crs2 = nn_dynamic_pointer_cast<ProjectedCRS>(obj2);
+    ASSERT_TRUE(crs2 != nullptr);
+
+    EXPECT_TRUE(crs->isEquivalentTo(
+        crs2.get(),
+        IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, projectedCRS_Krovak_EPSG_5221_as_PROJ_string) {
     auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
     auto crs = factory->createProjectedCRS("5221");
@@ -4750,6 +4883,18 @@ static DerivedGeographicCRSNNPtr createDerivedGeographicCRS() {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, derivedGeographicCRS_basic) {
+
+    auto derivedCRS = createDerivedGeographicCRS();
+    EXPECT_TRUE(derivedCRS->isEquivalentTo(derivedCRS.get()));
+    EXPECT_FALSE(derivedCRS->isEquivalentTo(
+        derivedCRS->baseCRS().get(), IComparable::Criterion::EQUIVALENT));
+    EXPECT_FALSE(derivedCRS->baseCRS()->isEquivalentTo(
+        derivedCRS.get(), IComparable::Criterion::EQUIVALENT));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, derivedGeographicCRS_WKT2) {
 
     auto expected = "GEODCRS[\"WMO Atlantic Pole\",\n"
@@ -4950,6 +5095,18 @@ static DerivedGeodeticCRSNNPtr createDerivedGeodeticCRS() {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, derivedGeodeticCRS_basic) {
+
+    auto derivedCRS = createDerivedGeodeticCRS();
+    EXPECT_TRUE(derivedCRS->isEquivalentTo(derivedCRS.get()));
+    EXPECT_FALSE(derivedCRS->isEquivalentTo(
+        derivedCRS->baseCRS().get(), IComparable::Criterion::EQUIVALENT));
+    EXPECT_FALSE(derivedCRS->baseCRS()->isEquivalentTo(
+        derivedCRS.get(), IComparable::Criterion::EQUIVALENT));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, derivedGeodeticCRS_WKT2) {
 
     auto expected = "GEODCRS[\"Derived geodetic CRS\",\n"
@@ -5796,6 +5953,31 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
                       CoordinateOperationContext::IntermediateCRSUse::NEVER),
                   crs_5340);
     }
+
+    // Check that we get the same result from an EPSG code and a CRS created
+    // from its WKT1 representation.
+    {
+        // Pulkovo 1942 / CS63 zone A2
+        auto crs = factory->createCoordinateReferenceSystem("2936");
+
+        // Two candidate transformations found, so not picking up any
+        EXPECT_EQ(crs->createBoundCRSToWGS84IfPossible(
+                      dbContext,
+                      CoordinateOperationContext::IntermediateCRSUse::NEVER),
+                  crs);
+
+        auto wkt = crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext)
+                .get());
+        auto obj =
+            WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
+        auto crs_from_wkt = nn_dynamic_pointer_cast<CRS>(obj);
+        ASSERT_TRUE(crs_from_wkt != nullptr);
+        EXPECT_EQ(crs_from_wkt->createBoundCRSToWGS84IfPossible(
+                      dbContext,
+                      CoordinateOperationContext::IntermediateCRSUse::NEVER),
+                  crs_from_wkt);
+    }
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp
index 1005d49b..ff86c4d3 100644
--- a/test/unit/test_factory.cpp
+++ b/test/unit/test_factory.cpp
@@ -531,8 +531,10 @@ TEST(factory, AuthorityFactory_createGeodeticCRS_geographic2D) {
     EXPECT_EQ(gcrs->identifiers()[0]->code(), "4326");
     EXPECT_EQ(*(gcrs->identifiers()[0]->codeSpace()), "EPSG");
     EXPECT_EQ(*(gcrs->name()->description()), "WGS 84");
-    EXPECT_TRUE(
-        gcrs->datum()->isEquivalentTo(factory->createDatum("6326").get()));
+    ASSERT_TRUE(gcrs->datum() == nullptr);
+    ASSERT_TRUE(gcrs->datumEnsemble() != nullptr);
+    EXPECT_TRUE(gcrs->datumEnsemble()->isEquivalentTo(
+        factory->createDatumEnsemble("6326").get()));
     EXPECT_TRUE(gcrs->coordinateSystem()->isEquivalentTo(
         factory->createCoordinateSystem("6422").get()));
     auto domain = crs->domains()[0];
@@ -566,8 +568,10 @@ TEST(factory, AuthorityFactory_createGeodeticCRS_geographic3D) {
     EXPECT_EQ(gcrs->identifiers()[0]->code(), "4979");
     EXPECT_EQ(*(gcrs->identifiers()[0]->codeSpace()), "EPSG");
     EXPECT_EQ(*(gcrs->name()->description()), "WGS 84");
-    EXPECT_TRUE(
-        gcrs->datum()->isEquivalentTo(factory->createDatum("6326").get()));
+    ASSERT_TRUE(gcrs->datum() == nullptr);
+    ASSERT_TRUE(gcrs->datumEnsemble() != nullptr);
+    EXPECT_TRUE(gcrs->datumEnsemble()->isEquivalentTo(
+        factory->createDatumEnsemble("6326").get()));
     EXPECT_TRUE(gcrs->coordinateSystem()->isEquivalentTo(
         factory->createCoordinateSystem("6423").get()));
 }
@@ -582,8 +586,10 @@ TEST(factory, AuthorityFactory_createGeodeticCRS_geocentric) {
     EXPECT_EQ(crs->identifiers()[0]->code(), "4978");
     EXPECT_EQ(*(crs->identifiers()[0]->codeSpace()), "EPSG");
     EXPECT_EQ(*(crs->name()->description()), "WGS 84");
-    EXPECT_TRUE(
-        crs->datum()->isEquivalentTo(factory->createDatum("6326").get()));
+    ASSERT_TRUE(crs->datum() == nullptr);
+    ASSERT_TRUE(crs->datumEnsemble() != nullptr);
+    EXPECT_TRUE(crs->datumEnsemble()->isEquivalentTo(
+        factory->createDatumEnsemble("6326").get()));
     EXPECT_TRUE(crs->coordinateSystem()->isEquivalentTo(
         factory->createCoordinateSystem("6500").get()));
 }
@@ -613,6 +619,20 @@ TEST(factory, AuthorityFactory_createVerticalCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(factory, AuthorityFactory_createVerticalCRS_with_datum_ensemble) {
+    auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    EXPECT_THROW(factory->createVerticalCRS("-1"),
+                 NoSuchAuthorityCodeException);
+
+    auto crs = factory->createVerticalCRS("9451"); // BI height
+    ASSERT_TRUE(crs->datum() == nullptr);
+    ASSERT_TRUE(crs->datumEnsemble() != nullptr);
+    EXPECT_TRUE(crs->datumEnsemble()->isEquivalentTo(
+        factory->createDatumEnsemble("1288").get()));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(factory, AuthorityFactory_createConversion) {
     auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
     EXPECT_THROW(factory->createConversion("-1"), NoSuchAuthorityCodeException);
@@ -3062,6 +3082,12 @@ TEST(factory, createObjectsFromName) {
                   .size(),
               1U);
 
+    // Exact name, but with other CRS that have an aliases to it ==> should
+    // match only the CRS with the given name, not those other CRS.
+    EXPECT_EQ(factory->createObjectsFromName("ETRS89 / UTM zone 32N", {}, false)
+                  .size(),
+              1U);
+
     // Prime meridian
     EXPECT_EQ(factoryEPSG->createObjectsFromName("Paris", {}, false, 2).size(),
               1U);
@@ -3166,6 +3192,29 @@ TEST(factory, createObjectsFromName) {
             .size(),
         1U);
 
+    {
+        auto res = factory->createObjectsFromName(
+            "World Geodetic System 1984 ensemble",
+            {AuthorityFactory::ObjectType::DATUM_ENSEMBLE}, false);
+        EXPECT_EQ(res.size(), 1U);
+        if (!res.empty()) {
+            EXPECT_EQ(res.front()->getEPSGCode(), 6326);
+            EXPECT_TRUE(dynamic_cast<DatumEnsemble *>(res.front().get()) !=
+                        nullptr);
+        }
+    }
+
+    {
+        auto res = factory->createObjectsFromName(
+            "World Geodetic System 1984 ensemble", {}, false);
+        EXPECT_EQ(res.size(), 1U);
+        if (!res.empty()) {
+            EXPECT_EQ(res.front()->getEPSGCode(), 6326);
+            EXPECT_TRUE(dynamic_cast<DatumEnsemble *>(res.front().get()) !=
+                        nullptr);
+        }
+    }
+
     const auto types = std::vector<AuthorityFactory::ObjectType>{
         AuthorityFactory::ObjectType::PRIME_MERIDIAN,
         AuthorityFactory::ObjectType::ELLIPSOID,
@@ -3187,6 +3236,7 @@ TEST(factory, createObjectsFromName) {
         AuthorityFactory::ObjectType::CONVERSION,
         AuthorityFactory::ObjectType::TRANSFORMATION,
         AuthorityFactory::ObjectType::CONCATENATED_OPERATION,
+        AuthorityFactory::ObjectType::DATUM_ENSEMBLE,
     };
     for (const auto type : types) {
         factory->createObjectsFromName("i_dont_exist", {type}, false, 1);
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 555d1159..4e903473 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -488,6 +488,62 @@ TEST(wkt_parse, wkt1_EPSG_4807_grad_mess) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt1_esri_EPSG_4901_grad) {
+    auto obj =
+        WKTParser()
+            .attachDatabaseContext(DatabaseContext::create())
+            .createFromWKT("GEOGCS[\"GCS_ATF_Paris\",DATUM[\"D_ATF\","
+                           "SPHEROID[\"Plessis_1817\",6376523.0,308.64]],"
+                           "PRIMEM[\"Paris_RGS\",2.33720833333333],"
+                           "UNIT[\"Grad\",0.0157079632679489]]");
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto datum = crs->datum();
+    auto primem = datum->primeMeridian();
+    EXPECT_EQ(primem->nameStr(), "Paris RGS");
+    // The PRIMEM is really in degree
+    EXPECT_EQ(primem->longitude().unit(), UnitOfMeasure::DEGREE);
+    EXPECT_NEAR(primem->longitude().value(), 2.33720833333333, 1e-14);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt2_epsg_org_EPSG_4901_PRIMEM_weird_sexagesimal_DMS) {
+    // Current epsg.org output may use the EPSG:9110 "sexagesimal DMS"
+    // unit and a DD.MMSSsss value, but this will likely be changed to
+    // use decimal degree.
+    auto obj = WKTParser().createFromWKT(
+        "GEOGCRS[\"ATF (Paris)\","
+        "  DATUM[\"Ancienne Triangulation Francaise (Paris)\","
+        "    ELLIPSOID[\"Plessis 1817\",6376523,308.64,"
+        "      LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]],"
+        "      ID[\"EPSG\",7027]],"
+        "    ID[\"EPSG\",6901]],"
+        "  PRIMEM[\"Paris RGS\",2.201395,"
+        "    ANGLEUNIT[\"sexagesimal DMS\",1,ID[\"EPSG\",9110]],"
+        "    ID[\"EPSG\",8914]],"
+        "  CS[ellipsoidal,2,"
+        "    ID[\"EPSG\",6403]],"
+        "  AXIS[\"Geodetic latitude (Lat)\",north,"
+        "    ORDER[1]],"
+        "  AXIS[\"Geodetic longitude (Lon)\",east,"
+        "    ORDER[2]],"
+        "  ANGLEUNIT[\"grad\",0.015707963267949,ID[\"EPSG\",9105]],"
+        "  USAGE[SCOPE[\"Geodesy.\"],AREA[\"France - mainland onshore.\"],"
+        "  BBOX[42.33,-4.87,51.14,8.23]],"
+        "ID[\"EPSG\",4901]]");
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto datum = crs->datum();
+    auto primem = datum->primeMeridian();
+    EXPECT_EQ(primem->longitude().unit(), UnitOfMeasure::DEGREE);
+    EXPECT_NEAR(primem->longitude().value(), 2.33720833333333, 1e-14);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, wkt1_geographic_old_datum_name_from_EPSG_code) {
     auto wkt =
         "GEOGCS[\"S-JTSK (Ferro)\",\n"
@@ -619,6 +675,52 @@ TEST(wkt_parse, wkt1_geographic_with_PROJ4_extension) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt1_geographic_epsg_org_api_4326) {
+    // Output from
+    // https://apps.epsg.org/api/v1/CoordRefSystem/4326/export/?format=wkt&formatVersion=1
+    // using a datum ensemble name
+    auto wkt =
+        "GEOGCS[\"WGS 84\",DATUM[\"World Geodetic System 1984 ensemble\","
+        "SPHEROID[\"WGS 84\",6378137,298.257223563,"
+        "AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],"
+        "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],"
+        "UNIT[\"degree (supplier to define representation)\","
+        "0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],"
+        "AXIS[\"Lat\",north],AXIS[\"Lon\",east],"
+        "AUTHORITY[\"EPSG\",\"4326\"]]";
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto datum = crs->datum();
+    EXPECT_EQ(datum->nameStr(), "World Geodetic System 1984");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_geographic_epsg_org_api_4258) {
+    // Output from
+    // https://apps.epsg.org/api/v1/CoordRefSystem/4258/export/?format=wkt&formatVersion=1
+    // using a datum ensemble name
+    auto wkt = "GEOGCS[\"ETRS89\","
+               "DATUM[\"European Terrestrial Reference System 1989 ensemble\","
+               "SPHEROID[\"GRS 1980\",6378137,298.257222101,"
+               "AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6258\"]],"
+               "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],"
+               "UNIT[\"degree (supplier to define representation)\","
+               "0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],"
+               "AXIS[\"Lat\",north],AXIS[\"Lon\",east],"
+               "AUTHORITY[\"EPSG\",\"4258\"]]";
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto datum = crs->datum();
+    EXPECT_EQ(datum->nameStr(), "European Terrestrial Reference System 1989");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, wkt1_geocentric_with_PROJ4_extension) {
     auto wkt = "GEOCCS[\"WGS 84\",\n"
                "    DATUM[\"unknown\",\n"
@@ -1658,16 +1760,21 @@ TEST(wkt_parse, wkt2_projected) {
                "                ID[\"EPSG\",9122]],\n"
                "            ID[\"EPSG\",8801]],\n"
                "        PARAMETER[\"Longitude of natural origin\",3,\n"
-               "            ANGLEUNIT[\"degree\",0.0174532925199433,\n"
-               "                ID[\"EPSG\",9122]],\n"
+               // Volontary omit LENGTHUNIT to check the WKT grammar accepts
+               // Check that we default to degree
+               //"            ANGLEUNIT[\"degree\",0.0174532925199433,\n"
+               //"                ID[\"EPSG\",9122]],\n"
                "            ID[\"EPSG\",8802]],\n"
                "        PARAMETER[\"Scale factor at natural origin\",0.9996,\n"
-               "            SCALEUNIT[\"unity\",1,\n"
-               "                ID[\"EPSG\",9201]],\n"
+               // Check that we default to unity
+               //"            SCALEUNIT[\"unity\",1,\n"
+               //"                ID[\"EPSG\",9201]],\n"
                "            ID[\"EPSG\",8805]],\n"
                "        PARAMETER[\"False easting\",500000,\n"
-               "            LENGTHUNIT[\"metre\",1,\n"
-               "                ID[\"EPSG\",9001]],\n"
+               // Volontary omit LENGTHUNIT to check the WKT grammar accepts
+               // Check that we default to metre
+               //"            LENGTHUNIT[\"metre\",1,\n"
+               //"                ID[\"EPSG\",9001]],\n"
                "            ID[\"EPSG\",8806]],\n"
                "        PARAMETER[\"False northing\",0,\n"
                "            LENGTHUNIT[\"metre\",1,\n"
@@ -1879,6 +1986,60 @@ TEST(wkt_parse, wkt2_2019_projected_utm_3D) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt2_2019_projected_with_base_geocentric) {
+    auto wkt =
+        "PROJCRS[\"EPSG topocentric example B\",\n"
+        "    BASEGEODCRS[\"WGS 84\",\n"
+        "        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n"
+        "            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G730)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G873)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n"
+        "            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                LENGTHUNIT[\"metre\",1]],\n"
+        "            ENSEMBLEACCURACY[2.0]],\n"
+        "        PRIMEM[\"Greenwich\",0,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "        ID[\"EPSG\",4978]],\n"
+        "    CONVERSION[\"EPSG topocentric example B\",\n"
+        "        METHOD[\"Geocentric/topocentric conversions\",\n"
+        "            ID[\"EPSG\",9836]],\n"
+        "        PARAMETER[\"Geocentric X of topocentric origin\",3771793.97,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8837]],\n"
+        "        PARAMETER[\"Geocentric Y of topocentric origin\",140253.34,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8838]],\n"
+        "        PARAMETER[\"Geocentric Z of topocentric origin\",5124304.35,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8839]]],\n"
+        "    CS[Cartesian,3],\n"
+        "        AXIS[\"topocentric East (U)\",east,\n"
+        "            ORDER[1],\n"
+        "            LENGTHUNIT[\"metre\",1]],\n"
+        "        AXIS[\"topocentric North (V)\",north,\n"
+        "            ORDER[2],\n"
+        "            LENGTHUNIT[\"metre\",1]],\n"
+        "        AXIS[\"topocentric height (W)\",up,\n"
+        "            ORDER[3],\n"
+        "            LENGTHUNIT[\"metre\",1]],\n"
+        "    USAGE[\n"
+        "        SCOPE[\"Example only (fictitious).\"],\n"
+        "        AREA[\"Description of the extent of the CRS.\"],\n"
+        "        BBOX[-90,-180,90,180]],\n"
+        "    ID[\"EPSG\",5820]]";
+    auto dbContext = DatabaseContext::create();
+    // Need a database so that EPSG:4978 is resolved
+    auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_TRUE(crs->baseCRS()->isGeocentric());
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, projected_angular_unit_from_primem) {
     auto obj = WKTParser().createFromWKT(
         "PROJCRS[\"NTF (Paris) / Lambert Nord France\",\n"
@@ -10444,6 +10605,14 @@ TEST(io, createFromUserInput) {
                  ParsingException);
     EXPECT_THROW(createFromUserInput("foobar + EGM96 height", dbContext),
                  ParsingException);
+
+    {
+        auto obj = createFromUserInput("World Geodetic System 1984 ensemble",
+                                       dbContext);
+        auto ensemble = nn_dynamic_pointer_cast<DatumEnsemble>(obj);
+        ASSERT_TRUE(ensemble != nullptr);
+        EXPECT_EQ(ensemble->identifiers().size(), 1U);
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -11139,6 +11308,155 @@ TEST(json_import, projected_crs) {
 
 // ---------------------------------------------------------------------------
 
+TEST(json_import, projected_crs_with_geocentric_base) {
+    auto json = "{\n"
+                "  \"$schema\": \"foo\",\n"
+                "  \"type\": \"ProjectedCRS\",\n"
+                "  \"name\": \"EPSG topocentric example B\",\n"
+                "  \"base_crs\": {\n"
+                "    \"name\": \"WGS 84\",\n"
+                "    \"datum_ensemble\": {\n"
+                "      \"name\": \"World Geodetic System 1984 ensemble\",\n"
+                "      \"members\": [\n"
+                "        {\n"
+                "          \"name\": \"World Geodetic System 1984 (Transit)\"\n"
+                "        },\n"
+                "        {\n"
+                "          \"name\": \"World Geodetic System 1984 (G730)\"\n"
+                "        },\n"
+                "        {\n"
+                "          \"name\": \"World Geodetic System 1984 (G873)\"\n"
+                "        },\n"
+                "        {\n"
+                "          \"name\": \"World Geodetic System 1984 (G1150)\"\n"
+                "        },\n"
+                "        {\n"
+                "          \"name\": \"World Geodetic System 1984 (G1674)\"\n"
+                "        },\n"
+                "        {\n"
+                "          \"name\": \"World Geodetic System 1984 (G1762)\"\n"
+                "        }\n"
+                "      ],\n"
+                "      \"ellipsoid\": {\n"
+                "        \"name\": \"WGS 84\",\n"
+                "        \"semi_major_axis\": 6378137,\n"
+                "        \"inverse_flattening\": 298.257223563\n"
+                "      },\n"
+                "      \"accuracy\": \"2.0\"\n"
+                "    },\n"
+                "    \"coordinate_system\": {\n"
+                "      \"subtype\": \"Cartesian\",\n"
+                "      \"axis\": [\n"
+                "        {\n"
+                "          \"name\": \"Geocentric X\",\n"
+                "          \"abbreviation\": \"X\",\n"
+                "          \"direction\": \"geocentricX\",\n"
+                "          \"unit\": \"metre\"\n"
+                "        },\n"
+                "        {\n"
+                "          \"name\": \"Geocentric Y\",\n"
+                "          \"abbreviation\": \"Y\",\n"
+                "          \"direction\": \"geocentricY\",\n"
+                "          \"unit\": \"metre\"\n"
+                "        },\n"
+                "        {\n"
+                "          \"name\": \"Geocentric Z\",\n"
+                "          \"abbreviation\": \"Z\",\n"
+                "          \"direction\": \"geocentricZ\",\n"
+                "          \"unit\": \"metre\"\n"
+                "        }\n"
+                "      ]\n"
+                "    },\n"
+                "    \"id\": {\n"
+                "      \"authority\": \"EPSG\",\n"
+                "      \"code\": 4978\n"
+                "    }\n"
+                "  },\n"
+                "  \"conversion\": {\n"
+                "    \"name\": \"EPSG topocentric example B\",\n"
+                "    \"method\": {\n"
+                "      \"name\": \"Geocentric/topocentric conversions\",\n"
+                "      \"id\": {\n"
+                "        \"authority\": \"EPSG\",\n"
+                "        \"code\": 9836\n"
+                "      }\n"
+                "    },\n"
+                "    \"parameters\": [\n"
+                "      {\n"
+                "        \"name\": \"Geocentric X of topocentric origin\",\n"
+                "        \"value\": 3771793.97,\n"
+                "        \"unit\": \"metre\",\n"
+                "        \"id\": {\n"
+                "          \"authority\": \"EPSG\",\n"
+                "          \"code\": 8837\n"
+                "        }\n"
+                "      },\n"
+                "      {\n"
+                "        \"name\": \"Geocentric Y of topocentric origin\",\n"
+                "        \"value\": 140253.34,\n"
+                "        \"unit\": \"metre\",\n"
+                "        \"id\": {\n"
+                "          \"authority\": \"EPSG\",\n"
+                "          \"code\": 8838\n"
+                "        }\n"
+                "      },\n"
+                "      {\n"
+                "        \"name\": \"Geocentric Z of topocentric origin\",\n"
+                "        \"value\": 5124304.35,\n"
+                "        \"unit\": \"metre\",\n"
+                "        \"id\": {\n"
+                "          \"authority\": \"EPSG\",\n"
+                "          \"code\": 8839\n"
+                "        }\n"
+                "      }\n"
+                "    ]\n"
+                "  },\n"
+                "  \"coordinate_system\": {\n"
+                "    \"subtype\": \"Cartesian\",\n"
+                "    \"axis\": [\n"
+                "      {\n"
+                "        \"name\": \"Topocentric East\",\n"
+                "        \"abbreviation\": \"U\",\n"
+                "        \"direction\": \"east\",\n"
+                "        \"unit\": \"metre\"\n"
+                "      },\n"
+                "      {\n"
+                "        \"name\": \"Topocentric North\",\n"
+                "        \"abbreviation\": \"V\",\n"
+                "        \"direction\": \"north\",\n"
+                "        \"unit\": \"metre\"\n"
+                "      },\n"
+                "      {\n"
+                "        \"name\": \"Topocentric height\",\n"
+                "        \"abbreviation\": \"W\",\n"
+                "        \"direction\": \"up\",\n"
+                "        \"unit\": \"metre\"\n"
+                "      }\n"
+                "    ]\n"
+                "  },\n"
+                "  \"scope\": \"Example only (fictitious).\",\n"
+                "  \"area\": \"Description of the extent of the CRS.\",\n"
+                "  \"bbox\": {\n"
+                "    \"south_latitude\": -90,\n"
+                "    \"west_longitude\": -180,\n"
+                "    \"north_latitude\": 90,\n"
+                "    \"east_longitude\": 180\n"
+                "  },\n"
+                "  \"id\": {\n"
+                "    \"authority\": \"EPSG\",\n"
+                "    \"code\": 5820\n"
+                "  }\n"
+                "}";
+    auto obj = createFromUserInput(json, nullptr);
+    auto pcrs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(pcrs != nullptr);
+    EXPECT_TRUE(pcrs->baseCRS()->isGeocentric());
+    EXPECT_EQ(pcrs->exportToJSON(&(JSONFormatter::create()->setSchema("foo"))),
+              json);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(json_import, compound_crs) {
     auto json = "{\n"
                 "  \"$schema\": \"foo\",\n"
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index ce4b866b..d0085004 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4867,7 +4867,7 @@ TEST(operation, vertCRS_to_geogCRS_context) {
                 "3855"), // EGM2008 height
             authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
             ctxt);
-        ASSERT_EQ(list.size(), 2U);
+        ASSERT_EQ(list.size(), 3U);
         EXPECT_EQ(
             list[1]->exportToPROJString(
                 PROJStringFormatter::create(
@@ -4889,7 +4889,7 @@ TEST(operation, vertCRS_to_geogCRS_context) {
                 "3855"), // EGM2008 height
             authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
             ctxt);
-        ASSERT_EQ(list.size(), 2U);
+        ASSERT_EQ(list.size(), 3U);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
             "+proj=pipeline "
@@ -4938,7 +4938,7 @@ TEST(operation, vertCRS_to_geogCRS_context) {
             authFactory->createCoordinateReferenceSystem("7839"),
             // NZGD2000
             authFactory->createCoordinateReferenceSystem("4959"), ctxt);
-        ASSERT_EQ(list.size(), 1U);
+        ASSERT_EQ(list.size(), 2U);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
             "+proj=pipeline "
@@ -4964,7 +4964,7 @@ TEST(operation, vertCRS_to_geogCRS_context) {
             authFactory->createCoordinateReferenceSystem("8357"),
             // ETRS89
             authFactory->createCoordinateReferenceSystem("4937"), ctxt);
-        ASSERT_EQ(list.size(), 1U);
+        ASSERT_EQ(list.size(), 2U);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
             "+proj=pipeline "
@@ -8734,6 +8734,14 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
                   "+multiplier=1 "
                   "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
                   "+step +proj=axisswap +order=2,1");
+        EXPECT_EQ(list[0]->remarks(),
+                  "For NAD83(NSRS2007) to NAVD88 height (1) (EPSG:9173): Uses "
+                  "Geoid09 hybrid model. Replaced by 2012 model (CT code 6326)."
+                  "\n"
+                  "For NAD83(NSRS2007) to WGS 84 (1) (EPSG:15931): "
+                  "Approximation at the +/- 1m level assuming that "
+                  "NAD83(NSRS2007) is equivalent to WGS 84 within the accuracy "
+                  "of the transformation.");
     }
 
     // NAD83 + NAVD88 height --> WGS 84
@@ -8904,7 +8912,7 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) {
                                                                ctxt);
     // The checked value is not that important, but in case this changes,
     // likely due to a EPSG upgrade, worth checking
-    EXPECT_EQ(listCompoundToGeog2D.size(), 141U);
+    EXPECT_EQ(listCompoundToGeog2D.size(), 142U);
 
     auto listGeog2DToCompound =
         CoordinateOperationFactory::create()->createOperations(dst, nnSrc,
@@ -10025,6 +10033,192 @@ TEST(operation, createOperation_ossfuzz_18587) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, derivedGeographicCRS_with_to_wgs84_to_geographicCRS) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=ob_tran +o_proj=latlon +lat_0=0 +lon_0=180 +o_lat_p=18.0 "
+        "+o_lon_p=-200.0 +ellps=WGS84 +towgs84=1,2,3 +type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=WGS84 +type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+
+    {
+        auto op = CoordinateOperationFactory::create()->createOperation(
+            NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+        ASSERT_TRUE(op != nullptr);
+        std::string pipeline(
+            "+proj=pipeline "
+            "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+            "+step +inv +proj=ob_tran +o_proj=latlon +lat_0=0 +lon_0=180 "
+            "+o_lat_p=18 +o_lon_p=-200 +ellps=WGS84 "
+            "+step +proj=push +v_3 "
+            "+step +proj=cart +ellps=WGS84 "
+            "+step +proj=helmert +x=1 +y=2 +z=3 "
+            "+step +inv +proj=cart +ellps=WGS84 "
+            "+step +proj=pop +v_3 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+        EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+                  pipeline);
+
+        auto op2 = CoordinateOperationFactory::create()->createOperation(
+            NN_CHECK_ASSERT(src),
+            nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4326));
+        ASSERT_TRUE(op2 != nullptr);
+        EXPECT_EQ(op2->exportToPROJString(PROJStringFormatter::create().get()),
+                  pipeline + " +step +proj=axisswap +order=2,1");
+    }
+
+    {
+        auto op = CoordinateOperationFactory::create()->createOperation(
+            NN_CHECK_ASSERT(dst), NN_CHECK_ASSERT(src));
+        ASSERT_TRUE(op != nullptr);
+        std::string pipeline(
+            "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+            "+step +proj=push +v_3 "
+            "+step +proj=cart +ellps=WGS84 "
+            "+step +proj=helmert +x=-1 +y=-2 +z=-3 "
+            "+step +inv +proj=cart +ellps=WGS84 "
+            "+step +proj=pop +v_3 "
+            "+step +proj=ob_tran +o_proj=latlon +lat_0=0 +lon_0=180 "
+            "+o_lat_p=18 +o_lon_p=-200 +ellps=WGS84 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+        EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+                  "+proj=pipeline " + pipeline);
+
+        auto op2 = CoordinateOperationFactory::create()->createOperation(
+            nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4326),
+            NN_CHECK_ASSERT(src));
+        ASSERT_TRUE(op2 != nullptr);
+        EXPECT_EQ(op2->exportToPROJString(PROJStringFormatter::create().get()),
+                  "+proj=pipeline +step +proj=axisswap +order=2,1 " + pipeline);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, geographic_topocentric) {
+    auto wkt =
+        "PROJCRS[\"EPSG topocentric example A\",\n"
+        "    BASEGEOGCRS[\"WGS 84\",\n"
+        "        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n"
+        "            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G730)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G873)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n"
+        "            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                LENGTHUNIT[\"metre\",1]],\n"
+        "            ENSEMBLEACCURACY[2.0]],\n"
+        "        PRIMEM[\"Greenwich\",0,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "        ID[\"EPSG\",4979]],\n"
+        "    CONVERSION[\"EPSG topocentric example A\",\n"
+        "        METHOD[\"Geographic/topocentric conversions\",\n"
+        "            ID[\"EPSG\",9837]],\n"
+        "        PARAMETER[\"Latitude of topocentric origin\",55,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "            ID[\"EPSG\",8834]],\n"
+        "        PARAMETER[\"Longitude of topocentric origin\",5,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "            ID[\"EPSG\",8835]],\n"
+        "        PARAMETER[\"Ellipsoidal height of topocentric origin\",0,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8836]]],\n"
+        "    CS[Cartesian,3],\n"
+        "        AXIS[\"topocentric East (U)\",east,\n"
+        "            ORDER[1],\n"
+        "            LENGTHUNIT[\"metre\",1]],\n"
+        "        AXIS[\"topocentric North (V)\",north,\n"
+        "            ORDER[2],\n"
+        "            LENGTHUNIT[\"metre\",1]],\n"
+        "        AXIS[\"topocentric height (W)\",up,\n"
+        "            ORDER[3],\n"
+        "            LENGTHUNIT[\"metre\",1]],\n"
+        "    USAGE[\n"
+        "        SCOPE[\"Example only (fictitious).\"],\n"
+        "        AREA[\"Description of the extent of the CRS.\"],\n"
+        "        BBOX[-90,-180,90,180]],\n"
+        "    ID[\"EPSG\",5819]]";
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto dst = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(dst != nullptr);
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        GeographicCRS::EPSG_4979, NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
+              "+step +proj=cart +ellps=WGS84 "
+              "+step +proj=topocentric +lat_0=55 +lon_0=5 +h_0=0 +ellps=WGS84");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, geocentric_topocentric) {
+    auto wkt =
+        "PROJCRS[\"EPSG topocentric example B\",\n"
+        "    BASEGEODCRS[\"WGS 84\",\n"
+        "        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n"
+        "            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G730)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G873)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n"
+        "            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n"
+        "            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                LENGTHUNIT[\"metre\",1]],\n"
+        "            ENSEMBLEACCURACY[2.0]],\n"
+        "        PRIMEM[\"Greenwich\",0,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "        ID[\"EPSG\",4978]],\n"
+        "    CONVERSION[\"EPSG topocentric example B\",\n"
+        "        METHOD[\"Geocentric/topocentric conversions\",\n"
+        "            ID[\"EPSG\",9836]],\n"
+        "        PARAMETER[\"Geocentric X of topocentric origin\",3771793.97,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8837]],\n"
+        "        PARAMETER[\"Geocentric Y of topocentric origin\",140253.34,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8838]],\n"
+        "        PARAMETER[\"Geocentric Z of topocentric origin\",5124304.35,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8839]]],\n"
+        "    CS[Cartesian,3],\n"
+        "        AXIS[\"topocentric East (U)\",east,\n"
+        "            ORDER[1],\n"
+        "            LENGTHUNIT[\"metre\",1]],\n"
+        "        AXIS[\"topocentric North (V)\",north,\n"
+        "            ORDER[2],\n"
+        "            LENGTHUNIT[\"metre\",1]],\n"
+        "        AXIS[\"topocentric height (W)\",up,\n"
+        "            ORDER[3],\n"
+        "            LENGTHUNIT[\"metre\",1]],\n"
+        "    USAGE[\n"
+        "        SCOPE[\"Example only (fictitious).\"],\n"
+        "        AREA[\"Description of the extent of the CRS.\"],\n"
+        "        BBOX[-90,-180,90,180]],\n"
+        "    ID[\"EPSG\",5820]]";
+    auto dbContext = DatabaseContext::create();
+    // Need a database so that EPSG:4978 is resolved
+    auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
+    auto dst = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(dst != nullptr);
+    auto f(NS_PROJ::io::WKTFormatter::create(
+        NS_PROJ::io::WKTFormatter::Convention::WKT2_2019));
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        GeodeticCRS::EPSG_4978, NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=topocentric +X_0=3771793.97 +Y_0=140253.34 "
+              "+Z_0=5124304.35 +ellps=WGS84");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, mercator_variant_A_to_variant_B) {
     auto projCRS = ProjectedCRS::create(
         PropertyMap(), GeographicCRS::EPSG_4326,
@@ -11096,7 +11290,7 @@ TEST(operation, normalizeForVisualization) {
             src,
             authFactory->createCoordinateReferenceSystem("4979"), // WGS 84 3D
             ctxt);
-        ASSERT_EQ(list.size(), 2U);
+        ASSERT_EQ(list.size(), 3U);
         auto op = list[1];
         auto opNormalized = op->normalizeForVisualization();
         EXPECT_FALSE(opNormalized->_isEquivalentTo(op.get()));