Recently Proj4 and GDAL have been upgraded to enable vertical datum transformations using grid shift files. In Canada, the common official vertical datums are known as CGVD28 and CGVD2013. Realizations (geoid models, or grid shift files) are provided by NRCAN for these (HT_2.0.byn and CGG2013n83.byn, respectively) and other vertical datums.
For use in Proj4/GDAL/libLAS, etc., the .byn files must be converted to .gtx files. You can do this with a tool I wrote, called byn2gtx. The files are placed in the Proj4 data directory which, on Linux, is usually /usr/share/proj (if installed by a package manager) or /usr/local/share/proj (if installed from source). In any case, you’ll see the existing .gtx files in the directory.
For Proj4 and GDAL to find the grid shift files, they have to be registered in the vertcs.override.csv file in […]/share/gdal. This is what the entry looks like for CGVD2013n83 (EPSG:6647) and CGVD28 (EPSG:5713):
6647,CGVD2013 height,1127,Canadian Geodetic Vertical Datum of 2013,9001,1,0,6499,9665,CGG2013n83.gtx
5713,CGVD28 height,5114,Canadian Geodetic Vertical Datum of 1928,9001,1,0,6499,9665,HT2_0.gtx
Notice that the file name of the grid shift file is the last field in each line.
Now, if you wish to transform, say, a digital elevation raster from one datum to another, you can use gdalwarp Actually, you can’t use gdalwarp for this… my bad:
gdalwarp -s_srs EPSG:2956+5713 -t_srs EPSG:2956+6647 in.tif out.tif
This would convert a tif from NAD83(CSRS) UTM12, with CGVD28 orthometric heights, to NAD83(CRSR) UTM12 with CGVD2013 heights. The “+” character creates a compound coordinate reference system from a vertical and horizontal CRS. (Note: there is an equivalent compound CRS — 6655 — in the compdcrs.csv file in GDAL’s data directory, but it is currently the only Canadian compound CRS.)
You can let PostGIS know about these datums too, but PostGIS’ integer SRIDs can’t handle the compound CRSes that GDAL can. The solution is simple: define a new SRID for a compound CRS. For example, if I want an SRID to represent NAD83(CRSR) UTM12N, CGVD2013, I create a new SRID, say, 99991 and add it to the spatial_ref_sys table using the following query:
INSERT into spatial_ref_sys
(srid, auth_name, auth_srid, proj4text, srtext)
(99991, 'epsg', 2956,
'+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+geoidgrids=CGG2013n83.gtx +vunits=m +no_defs',
PROJCS["NAD83(CSRS) / UTM zone 12N",
VERT_DATUM["Canadian Geodetic Vertical Datum of 2013",2005,
Here, the first field is the SRID, 99991; the second is the authority, ‘epsg’; the third is the horizontal CRS’ SRID, 2956. The last two fields are the Proj4 text and WKT projection definition (these can be found on a registry such as epsg.io). You can see that the proj4text entry has a +geoidgrids argument that refers to the grid shift file, CGG2013n83.gtx. The VERT_CS element is also present in the WKT projection information.
Say you have some points in your PostGIS table with the SRID 2956 and you’ve configured another SRID, 99992, for UTM12N NAD83, CGVD28. To transform your dataset, you first apply the SRID for the compound CRS, then apply the transformation using the SQL below.
Converting from 2956 to 99991 does not imply a transformation, it merely relabels the geometry with the compound CRS. Converting from 99991 to 99992 does imply a transformation because the geoid references are different.
SELECT ST_Transform(ST_SetSRID(geom, 99991), 99992) FROM my_table;