[Proj] *SOLVED*: Creating a nadgrid shift file with GDAL

[Andre Joost]

Am 21.01.2013 15:12, schrieb Jan Hartmann:
> This is the solution, thanks Andre.

Glad I could help!

> To summarize: creating an ntv2 file with gdal utilities:
>
> Let's say you have a list of coordinates in a non-wgs84 projection, and
> a list of the same locations in wgs84. You want to create an ntv2-file
> to shift from one to the other.
>
> 1) transfer the original coordinates to latlon, without a wgs84 shift

Why that? If you want to use the self-defined CRS in Qgis, it should 
better be in "real" WGS84 coordinates, that is with the +towgs84 
parameters. After all, its the position of the bessel ellipsoid relative 
to the WGS84 ellipsoid we now all reference to. Your old CRS may have 
the same position, though.

> 2) compute for each point the distance in degrees between the wgs84
> location and the original latlon-location
... or take the values our ancestors have already computed long time 
ago, or the raster printed on the map. If the data is in meters, 
transform them like normal on the same ellipsoid.

> 3) separate the distance in its lon and lat component, and multiply both
> numbers by 3600
> 4) invert the lon-value (<==== this was the problem)

We should better say: exchange the sign. Not to conflict wih 1/x 
computation.

> 5) create two point files, one for the lons and one for the lats. The
> coordinates should be the wgs84 values, and the z-value the lon or lat.
> Look on the gdal_grid page how to create a vrt-version, it's also
> possible with postgis or any other ogr-format. Create also a point-file
> with the same coordinates and all z-values set to zero
> 6) create for these three point-files each a gridded raster with
> gdal_grid. Set reasonable boundaries in wgs84 with -txe and -ty, and
> choose output format Float32

Make sure all your data is within the boundaries, otherwise no shift is 
applied. ntv2 does no extrapolation, not even for a millimeter offset...

> 7) merge the grids into a four-band geotif: gdal_merge -separate -o
> <shiftgrid> <latgrid> <longrid> <zerogrid> <zerogrid>. Mind the order:
> lat first
> 8) invert the y-direction in <shiftgrid>: gdalwarp <shiftgrid>
> <new_shiftgrid> (<=== this was another problem. Is this a bug?)

I think there should be any kind of operation? Otherwise it should be 
done from the start. Just think of West being positive.


> 9) convert to ntv2: gdal_translate <new_shiftgrid> <ntv2_grid>
>
> I would appreciate if someone would test the procedure on their own
> data,

I will do in the near future with old maps from Germany.

> and if it works, to add it to the documentation.
Is there any, and where?


> PS I fully agree with Gerald Evenden on the reversal of the lon-values
> in the ntv2-format
>
> http://osgeo-org.1560.n6.nabble.com/NTv2-info-td3842825.html
>
> I seem to get the idea that the Canadians developed this format---comments
> please---and if they did then the designers deserve to be flogged for 8
> days
> with wet noodles for "positive east longitude." I can't believe it!!!! What
> is unbelievable arrogance by whoever had that stupid---yes, *VERY*
> stupid---idea.
>

Well, they did not forsee to have it be used world wide...

Greetings,
André Joost