[Proj] Dutch correction grid

[molnar at sas.elte.hu]

Dear Jan,

> I see. So this would be a floating point grid with two bands, the first
> for the x (lon) offset, the second for the y (lat) offset. All
> measurements in arcseconds.
>

I can share my experiences creating a grid shift binary (gsb) file for
Hungary (between "Hungarian Datum 1972 (HD72)" and WGS84). As NTv2
standard was invented on the Western Hemisphere, some parameters might
seem strange for those, who learned geodesy east from Greenwich...

0: If you have a list of coordinates of the identical points in your
starting and target geodetic systems in ellipsoid coordinates, you should
create two grids from this pointcloud.
(I used surfer for this) First grid should contain the latitude shifts,
and the second should contain the longitude shifts. I save my results in
.dat format. (Use an interpolation method, which calculates grid
node-point values outside the convex hull of your listed coordinates, or
add some virtual points near the corners of the created grid to ensure,
that you should have a valid transformation outside the web of your
control points, but inside the border your country.)
1: You should create a four column table. First and second columns should
be the latitude and longitude shift on the first datum grid points,
starting from the bottom right (!) corner of the grid in westward
direction. (I used excel to sort the list of points in the .dat files)
Third and fourth column should be the standard deviation of the latitude
and longitude shift value. (I think these values are never ever used in
the calculations. I used dummy values like 0.001000)
2: As longitude grows westward for the standard, longitude shifts should
be opposite sign as calculated (not lon(2nd datum)-lon(1st datum) but the
opposite).
3: the delimitation character for this table should not be 'tab' but
'white space'. (At least it worked for me...)
4: you should add a header and a one line footer to your table in a text
editor. The final ascii file should look like this: (see my comments after
# signs!)

NUM_OREC 11        # number of lines for the general header
NUM_SREC 11        # number of lines for the sub-header
NUM_FILE  1        #I had only one dataset (so only 1 sub-header is expected)
GS_TYPE SECONDS    #all angels are in arcseconds
VERSION NTv2.0     #
DATUM_F HD72        #name of your local datum, 1st datum
DATUM_T WGS84       #name of your 2nd datum
MAJOR_F 6378160.000  #1st datum ellipsoid half major axis in meter
MINOR_F 6356774.516  #1st ... minor axis
MAJOR_T 6378137.000  #2nd ... major ...
MINOR_T 6356752.314  #2nd ... minor ...
SUB_NAMEGGPSH95      #sub header begins: name of grid shift sub-dataset
PARENT  NONE         # ??? (do not change)
CREATED 04-10-10     #
UPDATED 04-10-10     #
S_LAT     164520.000000  # southernmost latitude in seconds
N_LAT     174960.000000  #northernmost ...
E_LONG    -82800.000000  #easternmost longitude in seconds. Mind the sign!
W_LONG    -57600.000000  #westernmost ... Mind the sign!
LAT_INC      360.000000  #latitude increment in seconds
LONG_INC     360.000000  #longitude increment in seconds
GS_COUNT 2130   #number of samples in the grid shift sub-dataset

# her come the data like:
-0.864000	 4.032000	0.001000	0.001000
-0.866400	 4.029600	0.001000	0.001000
-0.868800	 4.027200	0.001000	0.001000
-0.871200	 4.024800	0.001000	0.001000
...
...
# end line:
END  .33E+33

5. Save the file as something.asc
6. Use GDAit for converting the ascii file above to and .gsb file.
Download GDAit, install it, start it, and File->Import Grid file choose
Files of type: NTv2 ASCII grid (*.asc). Make this file as the default grid
shift file.
6. Create a .gsb grid using GDAit: File->Export Grid file and choose the
.gsb type.
7. Now you can copy and use this .gsb file for conversion in gdal
environment with nadgrids=... option.

> It's not meant for a datum shift, I've got two targets:
>
> - maps in epgs:28992 (Dutch coordinates, with corrections on a regular
> grid  up to 25 cm)
>
> - maps in the 19th century Bonne-based projection, with corrections for
> about 1000 triangulation points (mostly church towers), up to 60 meters.
> These should have to be gridded.

I suggest to create two separate .gsb files. One between the epsg:28992
system and WGS84 (or ETRS89) and one another file between the old dutch
grid and the WGS84 (or ETRS89) system.

Good luck!

Gabor Molnar