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

[Hermann Peifer]

On 2013-01-21 12:12, Jan Hartmann wrote:
>
> 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
> 2) compute for each point the distance in degrees between the wgs84
> location and the original latlon-location
> 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)
> 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
> 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?)
> 9) convert to ntv2: gdal_translate <new_shiftgrid> <ntv2_grid>
>
> I would appreciate if someone would test the procedure on their own
> data, and if it works, to add it to the documentation.
>
> Thanks for all your help!
>

I just successfully tested your procedure using the 31 test points given 
in BETA2007testdaten.csv [1].

And indeed: it is somewhat puzzling that one has to do step 8. I 
"solved" the issue by defining the grid's Y extents as:

gdal_grid (...) -tye ymax ymin # rather than providing the values in the 
documented order: -tye ymin ymax

Wild guess: This looks to me like a bug in gdal_grid.cpp, perhaps around 
line 1000.

Hermann

[1] http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/BETA2007testdaten.csv