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

[Peter Gagliardi]

Hi Even,

Thank you, your answer clears up a key confusion I had about the types of
grid shift files. It looks like I really meant .gtx wherever I said .gsb.
In my actual use case, I have files in a separate file format (.byn), which
I now assume will need to be converted to .gtx format to use with PROJ.

> You can use a BoundCRS for that.

Ah, I wasn’t aware of BOUNDCRS until now. To check my understanding, a
BoundCRS is one where you input coordinates in SOURCECRS (x, y, z). Then
when PROJ reads this in, it uses the ABRIDGEDTRANSFORMATION to convert to
TARGETCRS (x’, y’, z’) which becomes the effective source CRS. This can
then be converted +to some other CRS as usual. Is this an accurate
description?

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

This command helped point me in the right direction. I experimented using
+geoidgrids with other PROJ strings, and used projinfo to examine the WKT.
In my case, it looks like I will use something along these lines to support
2D EPSG + custom geoid file:

COMPOUNDCRS[“User-defined CRS”,

    <WKT for 2D horizontal CRS>

    BOUNDCRS[

         SOURCECRS[

            VERTCRS["unknown",

                VDATUM["unknown"],

                CS[vertical,1],

                    AXIS["gravity-related height (H)",up,

                        LENGTHUNIT["metre",1,

                            ID["EPSG",9001]]]]],
        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,3],

                    AXIS["latitude",north,

                        ORDER[1],

                        ANGLEUNIT["degree",0.0174532925199433]],

                    AXIS["longitude",east,

                        ORDER[2],

                        ANGLEUNIT["degree",0.0174532925199433]],

                    AXIS["ellipsoidal height",up,

                        ORDER[3],

                        LENGTHUNIT["metre",1]],

                ID["EPSG",4979]]],

   ABRIDGEDTRANSFORMATION["unknown to WGS84 ellipsoidal height",

            METHOD["GravityRelatedHeight to Geographic3D"],

            PARAMETERFILE["Geoid (height correction) model file","<GTX
filename>",

                ID["EPSG",8666]]]]]

I also have a related follow-up question about some of the PROJ parameters.
I now see that +geoidgrids=foo.gtx is for vertical shifts while
+nadgrids=foo.gsb is for horizontal shifts. If this is the case, when would
you use +proj=vgridshift +grids=foo.gtx?

Thank you,

Peter

On Fri, Dec 13, 2019 at 8:17 AM Even Rouault <even.rouault at spatialys.com>
wrote:

> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20191216/bfbc79bf/attachment-0001.html>