Horisontal and vertical gridshift drivers

[kbevers]

Until now gridshifts has not been working with the new API in proj.h since parsing of +nadgrids and +geoidgrids is build into pj_transform(). This PR introduces the possibility to do horizontal and vertical gridshifting with the pipeline API.

The vgridshift driver is a simple wrapper for pj_apply_vgridshift() and hgridshift is a wrapper for the functionality in pj_apply_gridshift.c. Horizontal gridshifting has historically been done by adding the +nadgrids parameter to the projection initializer string. With the changes in this PR it can now also be done in a projection pipeline by adding +proj=hgridshift +grids=grid1.gdb,grid2.gsb,gridN.gsb. Similarly vertical gridshifts can be done by using +proj=vgridshift +grids=grid1.gtx,grid2.gtx,gridN.gtx. With the old API you would use the +geoidgrids parameter.

I have deliberately used very generic names for the drivers since it can be expected to be used for more than just geoids and the North American Datum. With the generic names it is more obvious that the drivers can be used in many different cases, for instance transforming the vertical component from ellipsoidal heights to a chart datum like LAT.

A few examples comparing the old gridshifts style with the new:

# Old style horizontal gridshift
$ echo 173 -45 0 | ./cs2cs -f %3.15f +proj=latlong +datum=WGS84 +nadgrids=nzgd2kgrid0005.gsb
173.000107914443305     -44.998379535535008 0.000000000000000

# New style horizontal gridshift
# Input given in radians since I am using the driver with `cs2cs` which it is not really meant to be
$ echo 3.01941960595019 -0.7853981633974483 0 | ./cs2cs -f %3.15f +proj=hgridshift +grids=nzgd2kgrid0005.gsb
173.000107914443305     -44.998379535535008 0.000000000000000

# Old style vertical gridshift
$ echo 173 -45 0 | ./cs2cs +proj=latlon +geoidgrids=egm96_15.gtx
173dE   45dS 4.975

# New style vertical gridshift. Radians used for same reason as above.
$ echo 3.01941960595019 -0.7853981633974483 0 | ./cs2cs /cs2cs +proj=vgridshift +grids=egm96_15.gx
173dE   45dS 4.975

[kbevers]

I am now happy with this. @busstoptaktik would you please have a look and see if it makes sense to you?

At the moment the forward operation of hgridshift and vgridshift is an addition. It's just convention and as long as it is documented properly it doesn't matter, but still it might be a good idea to have a discussion of what is the most intuitive convention. For instance the vgridshift will probably mostly be used when going from ellipsoidal heights to orthometric heights, which is a subtraction of the geoid grid from the ellipsoidal heights. Taking that into account it could be more sensible to have vgridshift subtract the grid in the forward operation. Thoughts?

[rouault]

For instance the vgridshift will probably mostly be used when going from ellipsoidal heights to orthometric heights, which is a subtraction of the geoid grid from the ellipsoidal heights

Actually, my main use case is the reverse ;-) when orthorectifying images with RPC, the transform (long, lat, elev) -> (x, y) expects elevations to be relative to the WGS84 elliposoid, so you will go from orthometric to ellipsoidal.

I think it might be a good idea to keep the same conventions in the old and new API.

[kbevers]

I think it might be a good idea to keep the same conventions in the old and new API.

I agree. And that's how I've set it up so far, as demonstrated with the cs2cs calls above. But actually the underlying pj_apply_[v]gridshif-functions are called with the reverse conventions. That is, the forward functions in the new gridshift code calls pj_apply_[v]gridshift with inverse=1, and the reverse functions uses inverse=0. I believe I am doing the correct thing to stay consistent, but it is confusing nonetheless.

[busstoptaktik]

I agree with @rouault that keeping the sign convention consistent is preferable. I also agree with @kbevers that the convention is confusing for the vgridshift case.

But it may become slightly less confusing if we assume that the convention has been established in pre-GNSS times, which really is not that long ago, historically speaking.

Prior to GNSS, you could hardly measure anything but orthometric heights ("elevations"), so when needing to go from geodetic/polar to geocentric/cartesian coordinates, you would need a geoid model and convert your elevations to geometric ("ellipsoidal") heights.

Hence, going from the natural height system to the synthetic/geometrical would be seen as the forward direction.

This is in fine accordance with the convention of map projections: Forward is when projecting FROM the "natural" latitude/longitude system, TO the "synthetic" planar, rectangular map coordinates.

So the way to remember the sign conventions in Proj.4 may be that when you go from a (somewhat) natural system to a more planar geometrical, it is a step forward.

Which makes sense if you're a cartographer, and probably somewhat less sense if you are a geophysicist :-)

[busstoptaktik]

it -> is

[busstoptaktik]

A comment makes more sense at the block level. Consider:

/* Failure most likely means that the grid is missing */
if (0==P)
    return 10

[busstoptaktik]

cf comment in hgridshift about block level comments

[busstoptaktik]

Looks good to me. I only have a few comments on comment structure, which you abviously may ignore if you disagree.

I believe this is a highly valuable contribution, much improving the utility of the pipeline functionality.

You point out a potential oops in your initial presentation of the work (the one about needing to input radians in cs2cs). I believe this should be corrected in cs2cs, if at all.

[kbevers]

I have updated the comments in the code, and otherwise taken note of your comments about conventions. Could be useful at a later stage when writing proper documentation for this. Thanks for the review.

Which makes sense if you're a cartographer, and probably somewhat less sense if you are a geophysicist :-)

I can pretend to be a cartographer for a moment :-)

You point out a potential oops in your initial presentation of the work (the one about needing to input radians in cs2cs). I believe this should be corrected in cs2cs, if at all.

So how about adding a -D flag to cs2cs, which tells cs2cs to not convert the input to radians? -D for degrees, in case that wasn't obvious.

[kbevers]

curl -O http://download.osgeo.org/proj/proj-datumgrid-1.6.zip
'curl' is not recognized as an internal or external command,
operable program or batch file.
Command exited with code 1

So now curl is not available on AppVeyor!? I wonder what has changed... worked fine a few days ago.

[kbevers]

So now curl is not available on AppVeyor!? I wonder what has changed... worked fine a few days ago.

According to appveyor/ci#1426 git was updated which caused curl to disappear from PATH. This is being fixed with the next AppVeyor update. A ticket has been posted here. The update will happen on 2017-03-26 apparently. I'll merge this once AppVeyor is up and running again and all tests pass.