[Proj] Creating a nadgrid shift file
Andre Joost
andre+joost at nurfuerspam.de
Sun Jan 20 10:28:17 PST 2013
Hi,
I think I found a solution, if I have understood your test data correct:
What you have and what you want:
source:
id | lon | lat | deg E | deg N | dms E | dms N
----+--------+---------+-------+--------+----------+----------+
ll | 5.3876 | 52.1561 | 19395 | 187762 | 5d23'15" | 52d9'22"
lr | 5.4022 | 52.1561 | 19448 | 187762 | 5d24'08" | 52d9'22"
ur | 5.4022 | 52.1651 | 19448 | 187794 | 5d24'08" | 52d9'54"
ul | 5.3876 | 52.1651 | 19395 | 187794 | 5d23'15" | 52d9'54"
destination:
id | lon | lat | deg E | deg N | dms E | dms N
----+--------+---------+-------+--------+----------+----------+
ll | 5.37299| 52.14711| 19343 | 187730 | 5d22'23" | 52d8'50"
lr | 5.41681| 52.14711| 19501 | 187730 | 5d25'1" | 52d8'50"
ur | 5.41682| 52.17409| 19501 | 187826 | 5d25'1" | 52d10'27"
ul | 5.37298| 52.17409| 19343 | 187826 | 5d22'23" | 52d10'27"
I expanded the data a bit:
S 52.13811°N =187697"N
N 52.18309°N =187859"N
E -5.43142°E =-19553"E
W -5.35838°E =-19290"E
with zero values for the border cells
My ASCII ntv2 file has the following header:
NUM_OREC 11
NUM_SREC 11
NUM_FILE 1
GS_TYPE SECONDS
VERSION
SYSTEM_F
SYSTEM_T
MAJOR_F 0.000
MINOR_F 0.000
MAJOR_T 0.000
MINOR_T 0.000
SUB_NAME
PARENT NONE
CREATED
UPDATED
S_LAT 187697.000000
N_LAT 187859.000000
E_LONG -19553.000000
W_LONG -19290.000000
LAT_INC 32.400000
LONG_INC 52.600000
GS_COUNT 36
dx and dy are swapped in the file, upper right first, and dx values
inverted:
-32.370100 -52.59180 0.000000 0.000000
-32.363600 52.602400 0.000000 0.000000
32.35040 -52.64480 0.000000 0.000000
32.35610 52.634210 0.000000 0.000000
(and a lot of null value lines between)
Converted with GDAy, because the old GDAit won't run on Windows 7
anymore :-(
For cs2cs, I made an input file with your WGS84 values, and used
Amersfoort lat/lon instead of the projected CRS.
Put the following in a batch file:
echo epsg4326 RDgrid.txt >out.txt
cs2cs +init=epsg:4326 +to +proj=longlat +ellps=bessel
+nadgrids=D:\Karten\gdal\ntv2\RD.gsb +no_defs RD-grid.txt >>out.txt
echo Amersf ll nowgs84 - grid >>out.txt
cs2cs +proj=longlat +ellps=bessel +towgs84=0,0,0,0,0,0,0 +no_defs +to
+proj=longlat +ellps=bessel +nadgrids=D:\Karten\gdal\ntv2\RD.gsb
+no_defs RD-grid.txt >>out.txt
echo Amersf ll - grid >>out.txt
cs2cs +proj=longlat +ellps=bessel
+towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725
+no_defs +to +proj=longlat +ellps=bessel
+nadgrids=D:\Karten\gdal\ntv2\RD.gsb +no_defs RD-grid.txt >>out.txt
echo grid - Amersf ll nowgs84 >>out.txt
cs2cs +proj=longlat +ellps=bessel +nadgrids=D:\Karten\gdal\ntv2\RD.gsb
+no_defs +to +proj=longlat +ellps=bessel +towgs84=0,0,0,0,0,0,0 +no_defs
RD-grid.txt >>out.txt
echo grid 2 Amersf ll >>out.txt
cs2cs +proj=longlat +ellps=bessel +nadgrids=D:\Karten\gdal\ntv2\RD.gsb
+no_defs +to +proj=longlat +ellps=bessel
+towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725
+no_defs RD-grid.txt >>out.txt
gives this result:
epsg4326 RDgrid.txt
5d23'15.36"E 52d9'21.96"N 0.000
5d24'7.92"E 52d9'21.96"N 0.000
5d24'7.92"E 52d9'54.36"N 0.000
5d23'15.36"E 52d9'54.36"N 0.000
Amersf ll nowgs84 - grid
5d23'15.36"E 52d9'24.046"N -698.431
5d24'7.92"E 52d9'24.046"N -698.431
5d24'7.92"E 52d9'56.446"N -698.420
5d23'15.36"E 52d9'56.446"N -698.420
Amersf ll - grid
5d23'13.793"E 52d9'18.402"N 43.347
5d24'6.345"E 52d9'18.402"N 43.346
5d24'6.345"E 52d9'50.798"N 43.331
5d23'13.793"E 52d9'50.799"N 43.331
grid - Amersf ll nowgs84
5d22'23.077"E 52d8'47.83"N 698.441
5d25'0.512"E 52d8'47.823"N 698.441
5d25'0.565"E 52d10'24.624"N 698.411
5d22'23.046"E 52d10'24.63"N 698.411
grid 2 Amersf ll
5d22'24.637"E 52d8'53.47"N -43.360
5d25'2.095"E 52d8'53.464"N -43.359
5d25'2.148"E 52d10'30.276"N -43.313
5d22'24.606"E 52d10'30.282"N -43.314
Compared to the start, the "grid 2 Amersfoort" is the way to get the
desired coordinates. Don't ask me why...
Maybe a case of "If you give me rubbish, I do nothing"-behaviour of proj.
HTH,
André Joost
More information about the Proj
mailing list