[PROJ] Defining a WKT CRS with a custom grid shift file

[Even Rouault]

Peter,

> I am working with custom CRS transformations using cs2cs. For typical EPSG
> codes, this works well, but I now want to support custom grid shift files
> that are not listed in proj.db. I currently am unsure how to do this with
> WKT strings.
> 
> I want to describe the following two types of transformations:
> 
> 
>    1.
> 
>    Transform from horizontal CRS EPSG:XXXX to EPSG:4979 using yyyyy.gsb to
>    transform heights
>    2.
> 
>    Transform from EPSG:4979 to EPSG:XXXX using yyyyy.gsb to transform
>    heights
> 
> 
> Where yyyyy.gsb is a custom grid shift file. Is there a way to do this
> using a COMPOUNDCRS like the following?
> 
> COMPOUNDCRS[“User-Defined Projection”,
> 
>      <WKT for horizontal CRS EPSG:XXXX>,
>     VERTCRS[<whatever parameters are needed to specify yyyyy.gsb>]
> 
> ]

If it is a .gsb file, then it is normally a horizontal shift, not a vertical 
one. You shouldn't use a compoundCRs for that.

You can use a BoundCRS for that. Easy way to bootstrap one is to run

$ projinfo "+proj=longlat +ellps=GRS80 +nadgrids=foo.gsb +type=crs" \
		-o WKT2_2019 -q

BOUNDCRS[
    SOURCECRS[
        GEOGCRS["unknown",
            DATUM["Unknown based on GRS80 ellipsoid",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",7019]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8901]],
            CS[ellipsoidal,2],
                AXIS["longitude",east,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]],
                AXIS["latitude",north,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["unknown to WGS84",
        METHOD["NTv2",
            ID["EPSG",9615]],
        PARAMETERFILE["Latitude and longitude difference file","foo.gsb",
            ID["EPSG",8656]]]]

But if you just need it for cs2cs purposes, just use the above PROJ.4 style 
string.

Otherwise if it is really a vertical shift (should rather be a .gtx file), use 
"+proj=longlat +ellps=GRS80 +geoidgrids=foo.gtx +type=crs"

Even


-- 
Spatialys - Geospatial professional services
http://www.spatialys.com