pj_release: Use correct release date for 6.1.0

[Backport 6.1] Update to EPSG 9.6.2

NEWS

@@ -14,7 +14,7 @@
 
  o Update internal version of googletest to v1.8.1 (#1361)
 
- o Database update: EPSG v9.6.1, IGNF v3.0.3, ESRI 10.7.0
+ o Database update: EPSG v9.6.2 (#1462), IGNF v3.0.3, ESRI 10.7.0
    and add operation_version column (#1368)
 
  o Add proj_normalize_for_visualization() that attempts to apply axis

data/Makefile.am

@@ -64,7 +64,7 @@ proj.db: $(DATAPATH)/sql/*.sql
 	 for x in $(SQL_ORDERED_LIST); do \
 		export SQL_EXPANDED_LIST="$${SQL_EXPANDED_LIST} $(DATAPATH)/$$x"; \
 	 done; \
-	 if test "x$(PROJ_DB_CACHE_DIR)" != "x" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5" ; then \
+	 if test "x$(PROJ_DB_CACHE_DIR)" != "x" -a -x "$$(command -v md5sum)" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5" ; then \
 	    cat $${SQL_EXPANDED_LIST} | md5sum | diff - "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5" > /dev/null \
 	    && (echo "Reusing cached proj.db"; cp "$(PROJ_DB_CACHE_DIR)/proj.db" proj.db); \
 	 fi; \
@@ -84,7 +84,7 @@ proj.db: $(DATAPATH)/sql/*.sql
 		$(RM) proj.db; \
 		exit 1; \
 	 else \
-	   if test "x$(PROJ_DB_CACHE_DIR)" != "x" ; then \
+	   if test "x$(PROJ_DB_CACHE_DIR)" != "x" -a -x "$$(command -v md5sum)" ; then \
 		mkdir -p "$(PROJ_DB_CACHE_DIR)"; \
 		cat $${SQL_EXPANDED_LIST} | md5sum > "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5"; \
 		cp proj.db "$(PROJ_DB_CACHE_DIR)"; \

data/sql/customizations.sql

@@ -19,6 +19,25 @@ DELETE FROM "supersession" WHERE superseded_table_name = 'grid_transformation' A
                                  replacement_code = '8885' AND
                                  source = 'EPSG';
 
+-- ('EPSG','7001','ETRS89 to NAP height (1)') lacks an interpolationCRS with Amersfoort / EPSG:4289
+-- See https://salsa.debian.org/debian-gis-team/proj-rdnap/blob/debian/2008-8/Use%20of%20RDTRANS2008%20and%20NAPTRANS2008.pdf
+-- "The naptrans2008 VDatum-grid is referenced to the Bessel-1841 ellipsoid"
+CREATE TABLE dummy(foo);
+CREATE TRIGGER check_grid_transformation_epsg_7001
+BEFORE INSERT ON dummy
+FOR EACH ROW BEGIN
+    SELECT RAISE(ABORT, 'grid_transformation EPSG:7001 entry is not ETRS89 to NAP height (1)')
+        WHERE NOT EXISTS(SELECT 1 FROM grid_transformation WHERE auth_name = 'EPSG' AND code = '7001' AND name = 'ETRS89 to NAP height (1)');
+    SELECT RAISE(ABORT, 'grid_transformation EPSG:7001 entry has already an interpolationCRS')
+        WHERE EXISTS(SELECT 1 FROM grid_transformation WHERE auth_name = 'EPSG' AND code = '7001' AND interpolation_crs_auth_name IS NOT NULL);
+END;
+INSERT INTO dummy DEFAULT VALUES;
+DROP TRIGGER check_grid_transformation_epsg_7001;
+DROP TABLE dummy;
+UPDATE grid_transformation SET interpolation_crs_auth_name = 'EPSG',
+                               interpolation_crs_code = '4289'
+                           WHERE auth_name = 'EPSG' AND code = '7001';
+
 -- Define the allowed authorities, and their precedence, when researching a
 -- coordinate operation
 

data/sql/grid_alternatives.sql

@@ -798,3 +798,58 @@ INSERT INTO grid_alternatives(original_grid_name,
                               0,
                               'proj-datumgrid-oceania',
                               NULL, NULL, NULL, NULL);
+
+-- Netherlands / RDNAP (non-free grids)
+
+INSERT INTO grid_alternatives(original_grid_name,
+                              proj_grid_name,
+                              proj_grid_format,
+                              proj_method,
+                              inverse_direction,
+                              package_name,
+                              url, direct_download, open_license, directory)
+                      VALUES ('naptrans2008.gtx',
+                              'naptrans2008.gtx',
+                              'GTX',
+                              'vgridshift',
+                              0,
+                              NULL, -- package name
+                              'https://salsa.debian.org/debian-gis-team/proj-rdnap/raw/upstream/2008/naptrans2008.gtx',
+                              1, -- direct download
+                              0, -- non-freely licensed. See https://salsa.debian.org/debian-gis-team/proj-rdnap/raw/master/debian/copyright
+                              NULL);
+
+INSERT INTO grid_alternatives(original_grid_name,
+                              proj_grid_name,
+                              proj_grid_format,
+                              proj_method,
+                              inverse_direction,
+                              package_name,
+                              url, direct_download, open_license, directory)
+                      VALUES ('rdtrans2008.gsb',
+                              'rdtrans2008.gsb',
+                              'NTv2',
+                              'hgridshift',
+                              0,
+                              NULL, -- package name
+                              'https://salsa.debian.org/debian-gis-team/proj-rdnap/raw/upstream/2008/rdtrans2008.gsb',
+                              1, -- direct download
+                              0, -- non-freely licensed. See https://salsa.debian.org/debian-gis-team/proj-rdnap/raw/master/debian/copyright
+                              NULL);
+
+-- Belgium
+
+INSERT INTO grid_alternatives(original_grid_name,
+                              proj_grid_name,
+                              proj_grid_format,
+                              proj_method,
+                              inverse_direction,
+                              package_name,
+                              url, direct_download, open_license, directory)
+                      VALUES ('bd72lb72_etrs89lb08.gsb',
+                              'bd72lb72_etrs89lb08.gsb',
+                              'NTv2',
+                              'hgridshift',
+                              0,
+                              'proj-datumgrid-europe',
+                              NULL, NULL, NULL, NULL);

data/sql/grid_transformation.sql

@@ -256,8 +256,8 @@ INSERT INTO "grid_transformation" VALUES('EPSG','8364','S-JTSK [JTSK03] to S-JTS
 INSERT INTO "grid_transformation" VALUES('EPSG','8369','BD72 to ETRS89 (3)',NULL,NULL,'EPSG','9615','NTv2','EPSG','4313','EPSG','4258','EPSG','1347',0.01,'EPSG','8656','Latitude and longitude difference file','bd72lb72_etrs89lb08.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'IGN-Bel 0.01m',0);
 INSERT INTO "grid_transformation" VALUES('EPSG','8371','RGF93 to NGF IGN69 height (2)',NULL,NULL,'EPSG','1073','Geographic3D to GravityRelatedHeight (IGN2009)','EPSG','4965','EPSG','5720','EPSG','1326',0.02,'EPSG','8666','Geoid (height correction) model file','RAF09.mnt',NULL,NULL,NULL,NULL,NULL,NULL,'IGN Fra 09',0);
 INSERT INTO "grid_transformation" VALUES('EPSG','8372','RGF93 to IGN78 Corsica height (2)',NULL,NULL,'EPSG','1073','Geographic3D to GravityRelatedHeight (IGN2009)','EPSG','4965','EPSG','5721','EPSG','1327',0.05,'EPSG','8666','Geoid (height correction) model file','RAC09.mnt',NULL,NULL,NULL,NULL,NULL,NULL,'IGN Fra Cor 09',0);
-INSERT INTO "grid_transformation" VALUES('EPSG','8444','GDA94 to GDA2020 (4)',NULL,NULL,'EPSG','9615','NTv2','EPSG','4283','EPSG','7844','EPSG','4169',0.05,'EPSG','8656','Latitude and longitude difference file','XMAS_C_V1.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'ICSM-Cxr CI Conf',0);
-INSERT INTO "grid_transformation" VALUES('EPSG','8445','GDA94 to GDA2020 (5)',NULL,NULL,'EPSG','9615','NTv2','EPSG','4283','EPSG','7844','EPSG','1069',0.05,'EPSG','8656','Latitude and longitude difference file','COCOS_C_V1.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'ICSM-Cki Conf',0);
+INSERT INTO "grid_transformation" VALUES('EPSG','8444','GDA94 to GDA2020 (4)',NULL,NULL,'EPSG','9615','NTv2','EPSG','4283','EPSG','7844','EPSG','4169',0.05,'EPSG','8656','Latitude and longitude difference file','GDA94_GDA2020_conformal_christmas_island.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'ICSM-Cxr Conf',0);
+INSERT INTO "grid_transformation" VALUES('EPSG','8445','GDA94 to GDA2020 (5)',NULL,NULL,'EPSG','9615','NTv2','EPSG','4283','EPSG','7844','EPSG','1069',0.05,'EPSG','8656','Latitude and longitude difference file','GDA94_GDA2020_conformal_cocos_island.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'ICSM-Cck Conf',0);
 INSERT INTO "grid_transformation" VALUES('EPSG','8446','GDA94 to GDA2020 (3)',NULL,NULL,'EPSG','9615','NTv2','EPSG','4283','EPSG','7844','EPSG','2575',0.05,'EPSG','8656','Latitude and longitude difference file','GDA94_GDA2020_conformal.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'ICSM-Aus NTv2 Conf',0);
 INSERT INTO "grid_transformation" VALUES('EPSG','8447','GDA94 to GDA2020 (2)',NULL,NULL,'EPSG','9615','NTv2','EPSG','4283','EPSG','7844','EPSG','2575',0.05,'EPSG','8656','Latitude and longitude difference file','GDA94_GDA2020_conformal_and_distortion.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'ICSM-Aus Conf and Dist',0);
 INSERT INTO "grid_transformation" VALUES('EPSG','8451','GDA2020 to AHD height (1)',NULL,NULL,'EPSG','1048','Geographic3D to GravityRelatedHeight (Ausgeoid v2)','EPSG','7843','EPSG','5711','EPSG','4493',0.03,'EPSG','8666','Geoid (height correction) model file','AUSGeoid2020_windows_binary.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'GA-Aus 2020',0);

data/sql/helmert_transformation.sql

@@ -1089,8 +1089,10 @@ INSERT INTO "helmert_transformation" VALUES('EPSG','8892','LKS94 to ETRS89 (1)',
 INSERT INTO "helmert_transformation" VALUES('EPSG','8893','SRB_ETRS89 to ETRS89 (1)',NULL,NULL,'EPSG','9603','Geocentric translations (geog2D domain)','EPSG','8685','EPSG','4258','EPSG','3534',0.0,0.0,0.0,0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'IOGP-Sbr',0);
 INSERT INTO "helmert_transformation" VALUES('EPSG','8894','SRB_ETRS89 to WGS 84 (1)',NULL,NULL,'EPSG','9603','Geocentric translations (geog2D domain)','EPSG','8685','EPSG','4326','EPSG','3534',1.0,0.0,0.0,0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'IOGP-Srb',0);
 INSERT INTO "helmert_transformation" VALUES('EPSG','8895','CHTRF95 to ETRS89 (1)',NULL,NULL,'EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4151','EPSG','4258','EPSG','1286',0.0,0.0,0.0,0.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'IOGP-Che',0);
-INSERT INTO "helmert_transformation" VALUES('EPSG','8913','CR05 to CR-SIRGAS (1)',NULL,NULL,'EPSG','9606','Position Vector transformation (geog2D domain)','EPSG','5364','EPSG','8906','EPSG','1074',0.5,-0.16959,0.35312,0.51846,'EPSG','9001',-0.03385,0.16325,-0.03446,'EPSG','9104',0.03693,'EPSG','9202',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'IGN-Cri 2014',0);
-INSERT INTO "helmert_transformation" VALUES('EPSG','8914','CR05 to WGS 84 (2)',NULL,NULL,'EPSG','9606','Position Vector transformation (geog2D domain)','EPSG','5364','EPSG','4326','EPSG','1074',1.0,-0.16959,0.35312,0.51846,'EPSG','9001',-0.03385,0.16325,-0.03446,'EPSG','9104',0.03693,'EPSG','9202',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'IGN-Cri 2014',0);
+INSERT INTO "helmert_transformation" VALUES('EPSG','8913','CR05 to CR-SIRGAS (1)',NULL,NULL,'EPSG','9606','Position Vector transformation (geog2D domain)','EPSG','5364','EPSG','8906','EPSG','1074',0.5,-0.16959,0.35312,0.51846,'EPSG','9001',-0.03385,0.16325,-0.03446,'EPSG','9104',0.03693,'EPSG','9202',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'IGN-Cri 2014',1);
+INSERT INTO "helmert_transformation" VALUES('EPSG','8914','CR05 to WGS 84 (2)',NULL,NULL,'EPSG','9606','Position Vector transformation (geog2D domain)','EPSG','5364','EPSG','4326','EPSG','1074',1.0,-0.16959,0.35312,0.51846,'EPSG','9001',-0.03385,0.16325,-0.03446,'EPSG','9104',0.03693,'EPSG','9202',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'IGN-Cri 2014',1);
+INSERT INTO "helmert_transformation" VALUES('EPSG','8968','CR05 to CR-SIRGAS (1)',NULL,NULL,'EPSG','9606','Position Vector transformation (geog2D domain)','EPSG','5365','EPSG','8907','EPSG','1074',0.5,-0.16959,0.35312,0.51846,'EPSG','9001',-0.03385,0.16325,-0.03446,'EPSG','9104',0.03693,'EPSG','9202',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'IGN-Cri 2014',0);
+INSERT INTO "helmert_transformation" VALUES('EPSG','8969','CR05 to WGS 84 (2)',NULL,NULL,'EPSG','9606','Position Vector transformation (geog2D domain)','EPSG','5365','EPSG','4326','EPSG','1074',1.0,-0.16959,0.35312,0.51846,'EPSG','9001',-0.03385,0.16325,-0.03446,'EPSG','9104',0.03693,'EPSG','9202',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'IGN-Cri 2014',0);
 INSERT INTO "helmert_transformation" VALUES('EPSG','10085','Trinidad 1903 to WGS 84 (2)',NULL,NULL,'EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4302','EPSG','4326','EPSG','1339',3.0,-61.0,285.2,471.6,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'EOG-Tto Trin',0);
 INSERT INTO "helmert_transformation" VALUES('EPSG','10086','JAD69 to WGS 72 (1)',NULL,NULL,'EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4242','EPSG','4322','EPSG','3342',15.0,48.0,208.0,382.0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'SD-Jam',0);
 INSERT INTO "helmert_transformation" VALUES('EPSG','10089','Aratu to WGS 84 (5)',NULL,NULL,'EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4208','EPSG','4326','EPSG','2962',7.0,-163.466,317.396,-147.538,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'EXM-Bra Santos',0);

data/sql/metadata.sql

@@ -1,2 +1,2 @@
-INSERT INTO "metadata" VALUES('EPSG.VERSION', 'v9.6.1');
-INSERT INTO "metadata" VALUES('EPSG.DATE', '2019-03-25');
+INSERT INTO "metadata" VALUES('EPSG.VERSION', 'v9.6.2');
+INSERT INTO "metadata" VALUES('EPSG.DATE', '2019-05-08');

docs/source/geodesic.rst

@@ -93,7 +93,7 @@ The routines also calculate several other quantities of interest
   On a flat surface, we have :math:`m_{12}=s_{12}`.
 * :math:`M_{12}` and :math:`M_{21}` are geodesic scales.  If two
   geodesics are parallel at point 1 and separated by a small distance
-  :\math`dt`, then they are separated by a distance :math:`M_{12}\,dt` at
+  :math:`dt`, then they are separated by a distance :math:`M_{12}\,dt` at
   point 2.  :math:`M_{21}` is defined similarly (with the geodesics
   being parallel to one another at point 2).  :math:`M_{12}` and
   :math:`M_{21}` are dimensionless quantities.  On a flat surface,
@@ -159,7 +159,7 @@ areas is based on :cite:`Danielsen1989`.  These improve on the work of
 :cite:`Vincenty1975` in the following respects:
 
 * The results are accurate to round-off for terrestrial ellipsoids (the
-  error in the distance is less then 15 nanometers, compared to 0.1 mm
+  error in the distance is less than 15 nanometers, compared to 0.1 mm
   for Vincenty).
 * The solution of the inverse problem is always found.  (Vincenty's
   method fails to converge for nearly antipodal points.)

include/proj/coordinateoperation.hpp

@@ -1476,7 +1476,8 @@ class PROJ_GCC_DLL Transformation : public SingleOperation {
     PROJ_DLL static TransformationNNPtr
     createGravityRelatedHeightToGeographic3D(
         const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
-        const crs::CRSNNPtr &targetCRSIn, const std::string &filename,
+        const crs::CRSNNPtr &targetCRSIn, const crs::CRSPtr &interpolationCRSIn,
+        const std::string &filename,
         const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies);
 
     PROJ_DLL static TransformationNNPtr createVERTCON(

include/proj/internal/internal.hpp

@@ -140,6 +140,9 @@ std::string toupper(const std::string &osStr);
 PROJ_FOR_TEST std::vector<std::string> split(const std::string &osStr,
                                              char separator);
 
+PROJ_FOR_TEST std::vector<std::string> split(const std::string &osStr,
+                                             const std::string &separator);
+
 bool ci_equal(const char *a, const char *b) noexcept;
 
 #ifdef SUPPORT_DELETED_FUNCTION

include/proj/io.hpp

@@ -427,7 +427,8 @@ class PROJ_GCC_DLL PROJStringFormatter {
     PROJ_INTERNAL void setOmitProjLongLatIfPossible(bool omit);
     PROJ_INTERNAL bool omitProjLongLatIfPossible() const;
 
-    PROJ_INTERNAL void setOmitZUnitConversion(bool omit);
+    PROJ_INTERNAL void pushOmitZUnitConversion();
+    PROJ_INTERNAL void popOmitZUnitConversion();
     PROJ_INTERNAL bool omitZUnitConversion() const;
 
     PROJ_INTERNAL void setDropEarlyBindingsTerms(bool drop);

scripts/build_db.py

@@ -134,7 +134,7 @@ def fill_conversion(proj_db_cursor):
     trigger_sql = """
 CREATE TRIGGER conversion_method_check_insert_trigger
 BEFORE INSERT ON conversion
-FOR EACH ROW BEGIN
+BEGIN
 """
 
     proj_db_cursor.execute("SELECT coord_op_code, coord_op_name, area_of_use_code, coord_op_method_code, coord_op_method_name, epsg_coordoperation.deprecated FROM epsg.epsg_coordoperation LEFT JOIN epsg.epsg_coordoperationmethod USING (coord_op_method_code) WHERE coord_op_type = 'conversion' AND coord_op_name NOT LIKE '%to DMSH'")
@@ -598,8 +598,165 @@ proj_db_cursor = proj_db_conn.cursor()
 proj_db_cursor.execute('PRAGMA foreign_keys = 1;')
 
 ingest_sqlite_dump(proj_db_cursor, os.path.join(sql_dir_name, 'proj_db_table_defs.sql'))
+
+# A bit messy, but to avoid churn in our existing .sql files, we temporarily
+# recreate the original conversion and helmert_transformation tables
+# instead of the view in the true database.
+
+proj_db_cursor.execute("""DROP VIEW conversion;""")
+proj_db_cursor.execute("""DROP TABLE conversion_table;""")
+proj_db_cursor.execute("""CREATE TABLE conversion(
+    auth_name TEXT NOT NULL CHECK (length(auth_name) >= 1),
+    code TEXT NOT NULL CHECK (length(code) >= 1),
+    name TEXT NOT NULL CHECK (length(name) >= 2),
+
+    description TEXT,
+    scope TEXT,
+
+    area_of_use_auth_name TEXT NOT NULL,
+    area_of_use_code TEXT NOT NULL,
+
+    method_auth_name TEXT CHECK (method_auth_name IS NULL OR length(method_auth_name) >= 1),
+    method_code TEXT CHECK (method_code IS NULL OR length(method_code) >= 1),
+    method_name TEXT,
+
+    param1_auth_name TEXT,
+    param1_code TEXT,
+    param1_name TEXT,
+    param1_value FLOAT,
+    param1_uom_auth_name TEXT,
+    param1_uom_code TEXT,
+
+    param2_auth_name TEXT,
+    param2_code TEXT,
+    param2_name TEXT,
+    param2_value FLOAT,
+    param2_uom_auth_name TEXT,
+    param2_uom_code TEXT,
+
+    param3_auth_name TEXT,
+    param3_code TEXT,
+    param3_name TEXT,
+    param3_value FLOAT,
+    param3_uom_auth_name TEXT,
+    param3_uom_code TEXT,
+
+    param4_auth_name TEXT,
+    param4_code TEXT,
+    param4_name TEXT,
+    param4_value FLOAT,
+    param4_uom_auth_name TEXT,
+    param4_uom_code TEXT,
+
+    param5_auth_name TEXT,
+    param5_code TEXT,
+    param5_name TEXT,
+    param5_value FLOAT,
+    param5_uom_auth_name TEXT,
+    param5_uom_code TEXT,
+
+    param6_auth_name TEXT,
+    param6_code TEXT,
+    param6_name TEXT,
+    param6_value FLOAT,
+    param6_uom_auth_name TEXT,
+    param6_uom_code TEXT,
+
+    param7_auth_name TEXT,
+    param7_code TEXT,
+    param7_name TEXT,
+    param7_value FLOAT,
+    param7_uom_auth_name TEXT,
+    param7_uom_code TEXT,
+
+    deprecated BOOLEAN NOT NULL CHECK (deprecated IN (0, 1)),
+
+    CONSTRAINT pk_conversion PRIMARY KEY (auth_name, code)
+);""")
+
+proj_db_cursor.execute("""DROP VIEW helmert_transformation;""")
+proj_db_cursor.execute("""DROP TABLE helmert_transformation_table;""")
+proj_db_cursor.execute("""CREATE TABLE helmert_transformation(
+    auth_name TEXT NOT NULL CHECK (length(auth_name) >= 1),
+    code TEXT NOT NULL CHECK (length(code) >= 1),
+    name TEXT NOT NULL CHECK (length(name) >= 2),
+
+    description TEXT,
+    scope TEXT,
+
+    method_auth_name TEXT NOT NULL CHECK (length(method_auth_name) >= 1),
+    method_code TEXT NOT NULL CHECK (length(method_code) >= 1),
+    method_name NOT NULL CHECK (length(method_name) >= 2),
+
+    source_crs_auth_name TEXT NOT NULL,
+    source_crs_code TEXT NOT NULL,
+    target_crs_auth_name TEXT NOT NULL,
+    target_crs_code TEXT NOT NULL,
+
+    area_of_use_auth_name TEXT NOT NULL,
+    area_of_use_code TEXT NOT NULL,
+
+    accuracy FLOAT CHECK (accuracy >= 0),
+
+    tx FLOAT NOT NULL,
+    ty FLOAT NOT NULL,
+    tz FLOAT NOT NULL,
+    translation_uom_auth_name TEXT NOT NULL,
+    translation_uom_code TEXT NOT NULL,
+    rx FLOAT,
+    ry FLOAT,
+    rz FLOAT,
+    rotation_uom_auth_name TEXT,
+    rotation_uom_code TEXT,
+    scale_difference FLOAT,
+    scale_difference_uom_auth_name TEXT,
+    scale_difference_uom_code TEXT,
+    rate_tx FLOAT,
+    rate_ty FLOAT,
+    rate_tz FLOAT,
+    rate_translation_uom_auth_name TEXT,
+    rate_translation_uom_code TEXT,
+    rate_rx FLOAT,
+    rate_ry FLOAT,
+    rate_rz FLOAT,
+    rate_rotation_uom_auth_name TEXT,
+    rate_rotation_uom_code TEXT,
+    rate_scale_difference FLOAT,
+    rate_scale_difference_uom_auth_name TEXT,
+    rate_scale_difference_uom_code TEXT,
+    epoch FLOAT,
+    epoch_uom_auth_name TEXT,
+    epoch_uom_code TEXT,
+    px FLOAT, -- Pivot / evaluation point for Molodensky-Badekas
+    py FLOAT,
+    pz FLOAT,
+    pivot_uom_auth_name TEXT,
+    pivot_uom_code TEXT,
+
+    operation_version TEXT, -- normally mandatory in OGC Topic 2 but optional here
+
+    deprecated BOOLEAN NOT NULL CHECK (deprecated IN (0, 1)),
+
+    CONSTRAINT pk_helmert_transformation PRIMARY KEY (auth_name, code)
+);""")
+
+proj_db_cursor.execute("SELECT name, sql FROM sqlite_master WHERE type = 'table' AND name = 'projected_crs'")
+for (name, sql) in proj_db_cursor.fetchall():
+    proj_db_cursor.execute("DROP TABLE " + name)
+    proj_db_cursor.execute(sql.replace('conversion_table', 'conversion'))
+
+proj_db_cursor.execute("SELECT name, sql FROM sqlite_master WHERE type = 'view'")
+for (name, sql) in proj_db_cursor.fetchall():
+    if 'conversion_table' in sql:
+        proj_db_cursor.execute("DROP VIEW " + name)
+        proj_db_cursor.execute(sql.replace('conversion_table', 'conversion'))
+    elif 'helmert_transformation_table' in sql:
+        proj_db_cursor.execute("DROP VIEW " + name)
+        proj_db_cursor.execute(sql.replace('helmert_transformation_table', 'helmert_transformation'))
+
 proj_db_cursor.execute("ATTACH DATABASE '%s' AS epsg;" % epsg_tmp_db_filename)
 
+
 fill_unit_of_measure(proj_db_cursor)
 fill_ellipsoid(proj_db_cursor)
 fill_area(proj_db_cursor)
@@ -647,7 +804,7 @@ for line in proj_db_conn.iterdump():
             f = open(os.path.join(sql_dir_name, table_name) + '.sql', 'wb')
             f.write("--- This file has been generated by scripts/build_db.py. DO NOT EDIT !\n\n".encode('UTF-8'))
             files[table_name] = f
-        f.write((line + '\n').encode('UTF-8'))
+        f.write((line + '\n').replace('BEFORE INSERT ON conversion', 'INSTEAD OF INSERT ON conversion').encode('UTF-8'))
 #f = files['coordinate_operation']
 #for row in non_imported_operations:
 #    f.write(("--- Non imported: " + str(row) + '\n').encode('UTF-8'))

scripts/reference_exported_symbols.txt

@@ -279,6 +279,7 @@ osgeo::proj::internal::ci_find(std::string const&, char const*)
 osgeo::proj::internal::c_locale_stod(std::string const&)
 osgeo::proj::internal::replaceAll(std::string const&, std::string const&, std::string const&)
 osgeo::proj::internal::split(std::string const&, char)
+osgeo::proj::internal::split(std::string const&, std::string const&)
 osgeo::proj::internal::tolower(std::string const&)
 osgeo::proj::internal::toString(double, int)
 osgeo::proj::io::AuthorityFactory::~AuthorityFactory()
@@ -618,7 +619,7 @@ osgeo::proj::operation::Transformation::createGeocentricTranslations(osgeo::proj
 osgeo::proj::operation::Transformation::createGeographic2DOffsets(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, std::vector<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> >, std::allocator<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> > > > const&)
 osgeo::proj::operation::Transformation::createGeographic2DWithHeightOffsets(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, std::vector<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> >, std::allocator<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> > > > const&)
 osgeo::proj::operation::Transformation::createGeographic3DOffsets(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, std::vector<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> >, std::allocator<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> > > > const&)
-osgeo::proj::operation::Transformation::createGravityRelatedHeightToGeographic3D(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, std::string const&, std::vector<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> >, std::allocator<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> > > > const&)
+osgeo::proj::operation::Transformation::createGravityRelatedHeightToGeographic3D(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, std::shared_ptr<osgeo::proj::crs::CRS> const&, std::string const&, std::vector<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> >, std::allocator<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> > > > const&)
 osgeo::proj::operation::Transformation::createLongitudeRotation(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, osgeo::proj::common::Angle const&)
 osgeo::proj::operation::Transformation::createMolodensky(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, double, double, double, double, double, std::vector<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> >, std::allocator<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> > > > const&)
 osgeo::proj::operation::Transformation::createNTv2(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::crs::CRS> > const&, std::string const&, std::vector<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> >, std::allocator<dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::metadata::PositionalAccuracy> > > > const&)

src/iso19111/c_api.cpp

@@ -30,10 +30,12 @@
 #define FROM_PROJ_CPP
 #endif
 
+#include <algorithm>
 #include <cassert>
 #include <cstdarg>
 #include <cstring>
 #include <map>
+#include <memory>
 #include <new>
 #include <utility>
 #include <vector>
@@ -151,6 +153,7 @@ static PJ *pj_obj_create(PJ_CONTEXT *ctx, const IdentifiedObjectNNPtr &objIn) {
     }
     auto pj = pj_new();
     if (pj) {
+        pj->ctx = ctx;
         pj->descr = "ISO-19111 object";
         pj->iso_obj = objIn;
     }
@@ -6793,6 +6796,59 @@ int proj_cs_get_axis_info(PJ_CONTEXT *ctx, const PJ *cs, int index,
  * nullptr in case of error
  */
 PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ *obj) {
+
+    SANITIZE_CTX(ctx);
+    if (!obj->alternativeCoordinateOperations.empty()) {
+        try {
+            auto pjNew = std::unique_ptr<PJ>(pj_new());
+            if (!pjNew)
+                return nullptr;
+            pjNew->ctx = ctx;
+            for (const auto &alt : obj->alternativeCoordinateOperations) {
+                auto co = dynamic_cast<const CoordinateOperation *>(
+                    alt.pj->iso_obj.get());
+                if (co) {
+                    double minxSrc = alt.minxSrc;
+                    double minySrc = alt.minySrc;
+                    double maxxSrc = alt.maxxSrc;
+                    double maxySrc = alt.maxySrc;
+                    double minxDst = alt.minxDst;
+                    double minyDst = alt.minyDst;
+                    double maxxDst = alt.maxxDst;
+                    double maxyDst = alt.maxyDst;
+
+                    auto l_sourceCRS = co->sourceCRS();
+                    auto l_targetCRS = co->targetCRS();
+                    if (l_sourceCRS && l_targetCRS) {
+                        const bool swapSource =
+                            l_sourceCRS
+                                ->mustAxisOrderBeSwitchedForVisualization();
+                        if (swapSource) {
+                            std::swap(minxSrc, minySrc);
+                            std::swap(maxxSrc, maxySrc);
+                        }
+                        const bool swapTarget =
+                            l_targetCRS
+                                ->mustAxisOrderBeSwitchedForVisualization();
+                        if (swapTarget) {
+                            std::swap(minxDst, minyDst);
+                            std::swap(maxxDst, maxyDst);
+                        }
+                    }
+                    pjNew->alternativeCoordinateOperations.emplace_back(
+                        minxSrc, minySrc, maxxSrc, maxySrc, minxDst, minyDst,
+                        maxxDst, maxyDst,
+                        pj_obj_create(ctx, co->normalizeForVisualization()),
+                        co->nameStr());
+                }
+            }
+            return pjNew.release();
+        } catch (const std::exception &e) {
+            proj_log_debug(ctx, __FUNCTION__, e.what());
+            return nullptr;
+        }
+    }
+
     auto co = dynamic_cast<const CoordinateOperation *>(obj->iso_obj.get());
     if (!co) {
         proj_log_error(ctx, __FUNCTION__, "Object is not a CoordinateOperation "

src/iso19111/coordinateoperation.cpp

@@ -6958,11 +6958,11 @@ TransformationNNPtr Transformation::createNTv2(
 static TransformationNNPtr _createGravityRelatedHeightToGeographic3D(
     const util::PropertyMap &properties, bool inverse,
     const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn,
-    const std::string &filename,
+    const crs::CRSPtr &interpolationCRSIn, const std::string &filename,
     const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
 
     return Transformation::create(
-        properties, sourceCRSIn, targetCRSIn, nullptr,
+        properties, sourceCRSIn, targetCRSIn, interpolationCRSIn,
         util::PropertyMap().set(
             common::IdentifiedObject::NAME_KEY,
             inverse ? INVERSE_OF + PROJ_WKT2_NAME_METHOD_HEIGHT_TO_GEOG3D
@@ -6981,17 +6981,20 @@ static TransformationNNPtr _createGravityRelatedHeightToGeographic3D(
  * At minimum the name should be defined.
  * @param sourceCRSIn Source CRS.
  * @param targetCRSIn Target CRS.
+ * @param interpolationCRSIn Interpolation CRS. (might be null)
  * @param filename GRID filename.
  * @param accuracies Vector of positional accuracy (might be empty).
  * @return new Transformation.
  */
 TransformationNNPtr Transformation::createGravityRelatedHeightToGeographic3D(
     const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
-    const crs::CRSNNPtr &targetCRSIn, const std::string &filename,
+    const crs::CRSNNPtr &targetCRSIn, const crs::CRSPtr &interpolationCRSIn,
+    const std::string &filename,
     const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
 
     return _createGravityRelatedHeightToGeographic3D(
-        properties, false, sourceCRSIn, targetCRSIn, filename, accuracies);
+        properties, false, sourceCRSIn, targetCRSIn, interpolationCRSIn,
+        filename, accuracies);
 }
 
 // ---------------------------------------------------------------------------
@@ -7302,8 +7305,20 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
     auto targetCRS = op->targetCRS();
     std::string name;
     if (!forwardName.empty()) {
-        if (starts_with(forwardName, INVERSE_OF)) {
-            name = forwardName.substr(INVERSE_OF.size());
+        if (starts_with(forwardName, INVERSE_OF) ||
+            forwardName.find(" + ") != std::string::npos) {
+            auto tokens = split(forwardName, " + ");
+            for (size_t i = tokens.size(); i > 0;) {
+                i--;
+                if (!name.empty()) {
+                    name += " + ";
+                }
+                if (starts_with(tokens[i], INVERSE_OF)) {
+                    name += tokens[i].substr(INVERSE_OF.size());
+                } else {
+                    name += INVERSE_OF + tokens[i];
+                }
+            }
         } else if (!sourceCRS || !targetCRS ||
                    forwardName != buildOpName(opType, sourceCRS, targetCRS)) {
             name = INVERSE_OF + forwardName;
@@ -8195,13 +8210,14 @@ TransformationNNPtr Transformation::substitutePROJAlternativeGridNames(
                 return createGravityRelatedHeightToGeographic3D(
                            createPropertiesForInverse(self.as_nullable().get(),
                                                       true, false),
-                           targetCRS(), sourceCRS(), projFilename,
-                           coordinateOperationAccuracies())
+                           targetCRS(), sourceCRS(), interpolationCRS(),
+                           projFilename, coordinateOperationAccuracies())
                     ->inverseAsTransformation();
             } else {
                 return createGravityRelatedHeightToGeographic3D(
                     createSimilarPropertiesTransformation(self), sourceCRS(),
-                    targetCRS(), projFilename, coordinateOperationAccuracies());
+                    targetCRS(), interpolationCRS(), projFilename,
+                    coordinateOperationAccuracies());
             }
         }
     }
@@ -10972,19 +10988,19 @@ struct MyPROJStringExportableHorizVertical final
         // cppcheck-suppress functionStatic
         _exportToPROJString(io::PROJStringFormatter *formatter) const override {
 
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
         horizTransform->_exportToPROJString(formatter);
 
         formatter->startInversion();
         geogDst->addAngularUnitConvertAndAxisSwap(formatter);
         formatter->stopInversion();
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
 
         verticalTransform->_exportToPROJString(formatter);
 
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
         geogDst->addAngularUnitConvertAndAxisSwap(formatter);
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
     }
 };
 
@@ -11016,7 +11032,7 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
         // cppcheck-suppress functionStatic
         _exportToPROJString(io::PROJStringFormatter *formatter) const override {
 
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
 
         opSrcCRSToGeogCRS->_exportToPROJString(formatter);
 
@@ -11024,17 +11040,17 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
         interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter);
         formatter->stopInversion();
 
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
 
         verticalTransform->_exportToPROJString(formatter);
 
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
 
         interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter);
 
         opGeogCRStoDstCRS->_exportToPROJString(formatter);
 
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
     }
 };
 
@@ -12135,6 +12151,37 @@ CoordinateOperationFactory::Private::createOperations(
             }
         }
 
+        auto vertCRSOfBaseOfBoundSrc =
+            dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
+        if (vertCRSOfBaseOfBoundSrc && hubSrcGeog &&
+            hubSrcGeog->coordinateSystem()->axisList().size() == 3 &&
+            geogDst->coordinateSystem()->axisList().size() == 3) {
+            auto opsFirst = createOperations(sourceCRS, hubSrc, context);
+            auto opsSecond = createOperations(hubSrc, targetCRS, context);
+            if (!opsFirst.empty() && !opsSecond.empty()) {
+                for (const auto &opFirst : opsFirst) {
+                    for (const auto &opLast : opsSecond) {
+                        // Exclude artificial transformations from the hub
+                        // to the target CRS
+                        if (!opLast->hasBallparkTransformation()) {
+                            try {
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {opFirst, opLast},
+                                            !allowEmptyIntersection));
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                            }
+                        }
+                    }
+                }
+                if (!res.empty()) {
+                    return res;
+                }
+            }
+        }
+
         return createOperations(boundSrc->baseCRS(), targetCRS, context);
     }
 
@@ -12349,7 +12396,8 @@ CoordinateOperationFactory::Private::createOperations(
         const auto &componentsSrc = compoundSrc->componentReferenceSystems();
         if (!componentsSrc.empty()) {
             std::vector<CoordinateOperationNNPtr> horizTransforms;
-            if (componentsSrc[0]->extractGeographicCRS()) {
+            auto srcGeogCRS = componentsSrc[0]->extractGeographicCRS();
+            if (srcGeogCRS) {
                 horizTransforms =
                     createOperations(componentsSrc[0], targetCRS, context);
             }
@@ -12363,11 +12411,61 @@ CoordinateOperationFactory::Private::createOperations(
                 for (const auto &horizTransform : horizTransforms) {
                     for (const auto &verticalTransform : verticalTransforms) {
 
-                        auto op = createHorizVerticalPROJBased(
-                            sourceCRS, targetCRS, horizTransform,
-                            verticalTransform);
+                        crs::GeographicCRSPtr interpolationGeogCRS;
+                        auto transformationVerticalTransform =
+                            dynamic_cast<const Transformation *>(
+                                verticalTransform.get());
+                        if (transformationVerticalTransform) {
+                            auto interpTransformCRS =
+                                transformationVerticalTransform
+                                    ->interpolationCRS();
+                            if (interpTransformCRS) {
+                                auto nn_interpTransformCRS =
+                                    NN_NO_CHECK(interpTransformCRS);
+                                if (dynamic_cast<const crs::GeographicCRS *>(
+                                        nn_interpTransformCRS.get())) {
+                                    interpolationGeogCRS =
+                                        util::nn_dynamic_pointer_cast<
+                                            crs::GeographicCRS>(
+                                            nn_interpTransformCRS);
+                                }
+                            }
+                        }
+                        bool done = false;
+                        if (interpolationGeogCRS &&
+                            (interpolationGeogCRS->_isEquivalentTo(
+                                 srcGeogCRS.get(),
+                                 util::IComparable::Criterion::EQUIVALENT) ||
+                             interpolationGeogCRS->is2DPartOf3D(
+                                 NN_NO_CHECK(srcGeogCRS.get())))) {
+                            auto srcToInterp = createOperations(
+                                componentsSrc[0],
+                                NN_NO_CHECK(interpolationGeogCRS), context);
+                            auto interpToCompoundHoriz = createOperations(
+                                NN_NO_CHECK(interpolationGeogCRS),
+                                componentsSrc[0], context);
+                            if (!srcToInterp.empty() &&
+                                !interpToCompoundHoriz.empty()) {
+                                auto op = createHorizVerticalHorizPROJBased(
+                                    sourceCRS, componentsSrc[0],
+                                    srcToInterp.front(), verticalTransform,
+                                    interpToCompoundHoriz.front(),
+                                    interpolationGeogCRS);
+                                done = true;
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {op, horizTransform},
+                                            !allowEmptyIntersection));
+                            }
+                        }
+                        if (!done) {
+                            auto op = createHorizVerticalPROJBased(
+                                sourceCRS, targetCRS, horizTransform,
+                                verticalTransform);
 
-                        res.emplace_back(op);
+                            res.emplace_back(op);
+                        }
                     }
                 }
                 return res;

src/iso19111/internal.cpp

@@ -298,6 +298,21 @@ std::vector<std::string> split(const std::string &str, char separator) {
 
 // ---------------------------------------------------------------------------
 
+std::vector<std::string> split(const std::string &str,
+                               const std::string &separator) {
+    std::vector<std::string> res;
+    size_t lastPos = 0;
+    size_t newPos = 0;
+    while ((newPos = str.find(separator, lastPos)) != std::string::npos) {
+        res.push_back(str.substr(lastPos, newPos - lastPos));
+        lastPos = newPos + separator.size();
+    }
+    res.push_back(str.substr(lastPos));
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
 #ifdef _WIN32
 
 // For some reason, sqlite3_snprintf() in the sqlite3 builds used on AppVeyor

src/iso19111/io.cpp

@@ -3818,7 +3818,7 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
                     Transformation::createGravityRelatedHeightToGeographic3D(
                         PropertyMap().set(IdentifiedObject::NAME_KEY,
                                           transformationName),
-                        crs, GeographicCRS::EPSG_4979,
+                        crs, GeographicCRS::EPSG_4979, nullptr,
                         stripQuotes(extensionChildren[1]),
                         std::vector<PositionalAccuracyNNPtr>());
                 return nn_static_pointer_cast<CRS>(BoundCRS::create(
@@ -4922,7 +4922,7 @@ struct PROJStringFormatter::Private {
     };
     std::vector<InversionStackElt> inversionStack_{InversionStackElt()};
     bool omitProjLongLatIfPossible_ = false;
-    bool omitZUnitConversion_ = false;
+    std::vector<bool> omitZUnitConversion_{false};
     DatabaseContextPtr dbContext_{};
     bool useApproxTMerc_ = false;
     bool addNoDefs_ = true;
@@ -5939,15 +5939,21 @@ bool PROJStringFormatter::omitProjLongLatIfPossible() const {
 
 // ---------------------------------------------------------------------------
 
-void PROJStringFormatter::setOmitZUnitConversion(bool omit) {
-    assert(d->omitZUnitConversion_ ^ omit);
-    d->omitZUnitConversion_ = omit;
+void PROJStringFormatter::pushOmitZUnitConversion() {
+    d->omitZUnitConversion_.push_back(true);
+}
+
+// ---------------------------------------------------------------------------
+
+void PROJStringFormatter::popOmitZUnitConversion() {
+    assert(d->omitZUnitConversion_.size() > 1);
+    d->omitZUnitConversion_.pop_back();
 }
 
 // ---------------------------------------------------------------------------
 
 bool PROJStringFormatter::omitZUnitConversion() const {
-    return d->omitZUnitConversion_;
+    return d->omitZUnitConversion_.back();
 }
 
 // ---------------------------------------------------------------------------
@@ -6995,7 +7001,7 @@ PROJStringParser::Private::buildBoundOrCompoundCRSIfNeeded(int iStep,
             Transformation::createGravityRelatedHeightToGeographic3D(
                 PropertyMap().set(IdentifiedObject::NAME_KEY,
                                   "unknown to WGS84 ellipsoidal height"),
-                crs, GeographicCRS::EPSG_4979, geoidgrids,
+                crs, GeographicCRS::EPSG_4979, nullptr, geoidgrids,
                 std::vector<PositionalAccuracyNNPtr>());
         auto boundvcrs =
             BoundCRS::create(vcrs, GeographicCRS::EPSG_4979, transformation);

src/release.cpp

@@ -11,7 +11,7 @@ char const pj_release[] =
     STR(PROJ_VERSION_MAJOR)"."
     STR(PROJ_VERSION_MINOR)"."
     STR(PROJ_VERSION_PATCH)", "
-    "September 1st, 2019";
+    "May 15th, 2019";
 
 const char *pj_get_release() {
     return pj_release;

src/strerrno.cpp

@@ -70,7 +70,7 @@ pj_err_list[] = {
     "argument not numerical or out of range",                          /* -58 */
     "inconsistent unit type between input and output",                 /* -59 */
     "arguments are mutually exclusive",                                /* -60 */
-    "generic error of unknow origin",                                  /* -61 */
+    "generic error of unknown origin",                                 /* -61 */
 
     /* When adding error messages, remember to update ID defines in
        projects.h, and transient_error array in pj_transform                  */

test/unit/test_c_api.cpp

@@ -3247,4 +3247,69 @@ TEST_F(CApi, proj_normalize_for_visualization) {
               "+step +proj=utm +zone=31 +ellps=WGS84");
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_normalize_for_visualization_with_alternatives) {
+
+    auto P = proj_create_crs_to_crs(m_ctxt, "EPSG:4326", "EPSG:3003", nullptr);
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto Pnormalized = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper keeper_Pnormalized(Pnormalized);
+    ASSERT_NE(Pnormalized, nullptr);
+
+    {
+        PJ_COORD c;
+        // Approximately Roma
+        c.lpz.lam = 12.5;
+        c.lpz.phi = 42;
+        c.lpz.z = 0;
+        c = proj_trans(Pnormalized, PJ_FWD, c);
+        EXPECT_NEAR(c.xy.x, 1789912.46264783037, 1e-8);
+        EXPECT_NEAR(c.xy.y, 4655716.25402576849, 1e-8);
+        auto projstr = proj_pj_info(Pnormalized).definition;
+        ASSERT_NE(projstr, nullptr);
+        EXPECT_EQ(std::string(projstr),
+                  "proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad "
+                  "step proj=push v_3 step proj=cart ellps=WGS84 "
+                  "step inv proj=helmert x=-104.1 y=-49.1 z=-9.9 rx=0.971 "
+                  "ry=-2.917 rz=0.714 s=-11.68 convention=position_vector "
+                  "step inv proj=cart ellps=intl step proj=pop v_3 "
+                  "step proj=tmerc lat_0=0 lon_0=9 k=0.9996 x_0=1500000 "
+                  "y_0=0 ellps=intl");
+    }
+
+    {
+        PJ_COORD c;
+        // Approximately Roma
+        c.xyz.x = 1789912.46264783037;
+        c.xyz.y = 4655716.25402576849;
+        c.xyz.z = 0;
+        c = proj_trans(Pnormalized, PJ_INV, c);
+        EXPECT_NEAR(c.lp.lam, 12.5, 1e-8);
+        EXPECT_NEAR(c.lp.phi, 42, 1e-8);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_normalize_for_visualization_with_alternatives_reverse) {
+
+    auto P = proj_create_crs_to_crs(m_ctxt, "EPSG:3003", "EPSG:4326", nullptr);
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto Pnormalized = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper keeper_Pnormalized(Pnormalized);
+    ASSERT_NE(Pnormalized, nullptr);
+
+    PJ_COORD c;
+    // Approximately Roma
+    c.xyz.x = 1789912.46264783037;
+    c.xyz.y = 4655716.25402576849;
+    c.xyz.z = 0;
+    c = proj_trans(Pnormalized, PJ_FWD, c);
+    EXPECT_NEAR(c.lp.lam, 12.5, 1e-8);
+    EXPECT_NEAR(c.lp.phi, 42, 1e-8);
+}
+
 } // namespace

test/unit/test_operation.cpp

@@ -5903,6 +5903,35 @@ TEST(operation, boundCRS_with_basecrs_with_extent_to_geogCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    // ETRS89 3D
+    auto src = authFactory->createCoordinateReferenceSystem("4937");
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 "
+        "+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel "
+        "+nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m "
+        "+type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        src, NN_NO_CHECK(dst), ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +inv +proj=vgridshift +grids=naptrans2008.gtx "
+              "+multiplier=1 "
+              "+step +inv +proj=hgridshift +grids=rdtrans2008.gsb "
+              "+step +proj=sterea +lat_0=52.1561605555556 "
+              "+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 "
+              "+y_0=463000 +ellps=bessel");
+}
+
+// ---------------------------------------------------------------------------
+
 static VerticalCRSNNPtr createVerticalCRS() {
     PropertyMap propertiesVDatum;
     propertiesVDatum.set(Identifier::CODESPACE_KEY, "EPSG")
@@ -5941,8 +5970,8 @@ static BoundCRSNNPtr createBoundVerticalCRS() {
     auto vertCRS = createVerticalCRS();
     auto transformation =
         Transformation::createGravityRelatedHeightToGeographic3D(
-            PropertyMap(), vertCRS, GeographicCRS::EPSG_4979, "egm08_25.gtx",
-            std::vector<PositionalAccuracyNNPtr>());
+            PropertyMap(), vertCRS, GeographicCRS::EPSG_4979, nullptr,
+            "egm08_25.gtx", std::vector<PositionalAccuracyNNPtr>());
     return BoundCRS::create(vertCRS, GeographicCRS::EPSG_4979, transformation);
 }
 
@@ -6727,6 +6756,22 @@ TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
     auto list =
         CoordinateOperationFactory::create()->createOperations(src, dst, ctxt);
     ASSERT_GE(list.size(), 1U);
+
+    {
+        // Important here is vgridshift before hgridshift
+        auto op_proj =
+            list[0]->exportToPROJString(PROJStringFormatter::create().get());
+        EXPECT_EQ(
+            op_proj,
+            "+proj=pipeline +step +inv +proj=sterea +lat_0=52.1561605555556 "
+            "+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 "
+            "+ellps=bessel "
+            "+step +proj=vgridshift +grids=naptrans2008.gtx +multiplier=1 "
+            "+step +proj=hgridshift +grids=rdtrans2008.gsb "
+            "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+            "+step +proj=axisswap +order=2,1");
+    }
+
     auto wkt2 =
         "COMPOUNDCRS[\"unknown\",\n"
         "  PROJCRS[\"unknown\",\n"
@@ -6759,8 +6804,12 @@ TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
     for (size_t i = 0; i < list.size(); i++) {
         const auto &op = list[i];
         const auto &op2 = list2[i];
-        EXPECT_TRUE(
-            op->isEquivalentTo(op2.get(), IComparable::Criterion::EQUIVALENT));
+        auto op_proj =
+            op->exportToPROJString(PROJStringFormatter::create().get());
+        auto op2_proj =
+            op2->exportToPROJString(PROJStringFormatter::create().get());
+        EXPECT_EQ(op_proj, op2_proj) << "op=" << op->nameStr()
+                                     << " op2=" << op2->nameStr();
     }
 }
 

travis/osx/before_install.sh

@@ -8,6 +8,7 @@ brew update
 brew install ccache
 brew install sqlite3
 brew install doxygen
+brew install md5sha1sum
 
 export PATH=$HOME/Library/Python/2.7/bin:$PATH
 # breathe=4.12.0 is the last version to work for us with sphinx 1.8.5 / Python 2