Update (July 2008):

Great news - the file provided on this page is no longer needed! Didier Richard from IGN has kindly provided an improved NTv2 version of the NTF-RGF93 gridshift file to the open source community and it is now available from the PROJ.4 website - download proj-datumgrid-1.4.zip and the new file ntf_r93.gsb is contained within that zip archive.

The new file is, as Didier explains, slightly more accurate than my file and provides centimetre precision. So I would recommend use of it rather than the file on this page.

The rest of the information in this page is now of historical interest only.

Creation of NTv2 Grid File for French NTF <--> RGF93 Datum Conversions

In June 2003 I commented on the GRASS mailing list about the inadequacy of the widely disseminated 3-parameter transformation from the traditional French NTF geodetic datum (Clark1880 ellipsoid) to the new ETRS89-consistent RGF93 datum (GRS1980 ellipsoid).

Date: Tue, 03 Jun 2003 00:32:39 +0200
From: Michel Wurtz 
To: Paul Kelly 
Subject: Re: [GRASS5] New Datums

Paul Kelly wrote:

> Another point is that it was very hard to find transformation parameters
> for France (I tried to add support for a few European countries but of
> course it is not complete). The above website only had one 3-parameter
> transformation for the whole of France which couldn't be very accurate.

Look at the official french survey site : all parameters and even the algorithms
used are described.  You can also load a grid af correction parameters to do
accurate transformation between NTF and GRS80 systems (same style as NAD to
NAD ?)

Look at http://www.ign.fr/ and more precisely at
   http://www.ign.fr/affiche_rubrique.asp?rbr_id=1068&lng_id=FR#44606
   http://www.ign.fr/affiche_rubrique.asp?rbr_id=1031&lng_id=FR
      
Everything is in french, but not too complicated as it is assentially values
and formulas.  Feel free to mail me directly if you have problems.

BTW, the systems used for overseas departements (French Guyana, Martinique, Guadeloupe,
Réunion are also described.  Docs are in html and formulas as pdf files.

You can also download Circé2000, a coordinate translator for Windows (alas) that will
give you a tool for verifying the Proj4 implementation and parameters...

-- 
Michel Wurtz - Auzeville-Tolosane

I thought it would be a useful project to convert the grid of projection parameters mentioned above into a format that could be used within the GRASS GIS for datum transformations. This meant transforming the file into a format compatible with the PROJ.4 system, which is used internally by GRASS to handle co-ordinate system conversions.

I decided I would try to create a file in the Canadian NTv2 format, as it is supported by PROJ.4 and I found many references to it on the internet. This suggested to me that it might be well on its way to becoming a de facto standard, and I thought I might as well do my bit to help it along.

Canadian file format

The best reference I found to this format was actually from the Department of Geomatics at the University of Melbourne in Australia (the NTv2 format is also officially used in Australia). They have produced some Windows software, GDAit, that can convert an NTv2 file between the ASCII and binary formats, and the detailed structure of the files is described quite adequately in appendix C of the GDAit User's Guide.

Link: University of Melbourne GDAit page

The main point to understand about the file is that it contains a grid of latitude and longitude co-ordinates, in the original (i.e. localised country) datum, and a latitude and longitude shift value to convert these lat/long co-ordinates into the co-ordinates (of the same real-world point) on the WGS84-compatible datum. This is fairly easy to understand intuitively. The latitude and longitude aren't given explicitly on each line of the file, only in the header. Each pair of shift values in the grid is listed row-wise and column-wise starting at the lower right corner of the grid and working backwards to the upper left. The values are given in seconds and on the same line a pair of accuracy / resolution figures for each shift value is also given. Therefore after the header there are four numbers on each line in the file. See the GDAit User's Guide above for more information.

French file format

The most useful webpage I have found regarding this format is the Institut Géographique National (IGN) download (téléchargement) page in France, in particular the User's Guide for the datum conversion grid file (Notice d'utilisation de la grille de paramètres de transformation de coordonnées GR3DF97A) and the file itself.

The French file format explicitly states the latitude and longitude on each line of the file. However these are lat/long co-ordinates in the WGS84-compatible datum, not in the local datum as in the Canadian file. Instead of giving latitude and longitude shifts, the French file gives the shifts in terms of cartesian geocentric co-ordinates, i.e. there are three shift values: dx, dy and dz. It is very important to note that these are given in the conventional `towgs84' shift direction. This means that to convert from the WGS84-compatible datum to the older NTF datum, the shift values must be SUBTRACTED from the WGS84/GRS80 geocentric co-ordinates. I overlooked this point (partly because of my incomplete understanding of the French instructions in the User's guide) and spent a lot of time wondering what was wrong later on.

The file also contains an accuracy figure for the transformation at each point on the grid. As far as I can understand this is given in metres on the ground, i.e. the projected distance in whatever zone is in use at that location

Transforming the Data

To convert the French ``backward dx dy dz shift for each point in the new datum'' to the Canadian ``forward latitude longitude shift for each point in the old datum'', the following procedure must be followed:

1. Transform grs80 latitude / longitude to grs80 geocentric x, y, z
2. SUBSTRACT the dx dy dz offset to give clark80 geocentric x, y, z
3. Transform clark80 geocentric to clark80 latitude / longitude and store the result
4. Subtract the clark80 latitude/longitude co-ordinates from the grs80 ones to give the clark80 to grs80 shift, and store these values along with the clark80 lat/long values

The above procedure can be accomplished on Unix using the cs2cs program from the PROJ.4 distribution and the awk text processing program. E.g. the first line in the French file is for latitude 41N, longitude 5W and the shift parameters are -165.027, -67.100 and +315.813. The following command pipeline will output a line in the format ``clark80 lat., clark80 long., clark80->grs80 lat. shift, clark80->grs80 long. shift'':

echo '-5.5 41' | \
cs2cs -f %.3f +proj=latlong +ellps=GRS80 +to +proj=geocent +ellps=GRS80 | \
awk '{printf "%.3f\t%.3f\t%.3f\n",$1-(-165.027),$2-(-67.100),$3-(315.813)}' | \
cs2cs -f %.9f +proj=geocent +ellps=clrk80 +to +proj=latlong +ellps=clrk80 | \
awk '{printf "%.9f\t%.9f\t%.9f\t%.9f\n", $2, $1, 41-$2, -5.5-$1}'

41.000036342    -5.499018175    -0.000036342    -0.000981825

The above command line was embedded in a highly inefficient Perl script to process the entire French grid file and output a file containing the NTF (Clark80) latitude and longitude values together with the shifts and accuracy (not described here), in a form suitable for input to GRASS using the s.in.ascii command.

I created a latitude/longitude location in GRASS and imported the file, then set the region to the following settings:

projection: 3 (Latitude-Longitude) zone: 0 datum: ntf ellipsoid: clark80 north: 52:03N south: 40:57N west: 5:33W east: 10:03E nsres: 0:06 ewres: 0:06 rows: 111 cols: 156

This meant that the centre of each grid cell was located at the 0.1 degree intervals used in the grid shift file. Then I interpolated the latitude shifts and longitude shifts separately into two raster layers using s.surf.rst (not too sure if it's properly valid for lat/long locations, and it gave many warnings about taking too long to find the points despite them being evenly distributed across the whole region) but it gave good smooth results in the end as it always does.

As the centre of each cell corresponded to the correct grid point, a simple 'r.to.sites -a' followed by 's.out.ascii -d' gave two text files containing the latitude and longitude shift values in a regular grid of NTF / Clark80 latitude / longitude co-ordinates.

The Canadian format requires the shift values in seconds, so the following calls to awk did the conversion:

cat latdiff.txt | awk '{printf "%.6f\n", $3*3600}' > latsecs.txt
cat longdiff.txt | awk '{printf "%.6f\n", $3*3600}' > longsecs.txt

I manually created a header for the ASCII NTv2 file (see hints in the GDAit manual):

NUM_OREC 11 NUM_SREC 11 NUM_FILE 1 GS_TYPE SECONDS VERSION PK SYSTEM_FNTF SYSTEM_TRGF93 MAJOR_F 6378249.145 MINOR_F 6356514.870 MAJOR_T 6378137.000 MINOR_T 6356752.314 SUB_NAMEFRANCE PARENT NONE CREATED 05062003 UPDATED 11062003 S_LAT 147600.000000 N_LAT 187200.000000 E_LONG -36000.000000 W_LONG 19800.000000 LAT_INC 360.000000 LONG_INC 360.000000 GS_COUNT 17316

Most important to note here and in the next stage is that the Canadian format gives longitudes as positive to the west, so the sign must be changed, and also that the co-ordinates are all given in seconds (degrees/3600).

And the following few commands create the final ASCII file. I didn't make use of the accuracy information in the French file yet; I think it would take a better knowledge of projections than I currently possess to convert the projected metres values into separate latitude and longitude errors. I have just given it as 1 second for every value in the file; I don't think PROJ.4 uses this information so it should be all right, certainly for use in GRASS.

paste revlatsecs.txt revlongsecs.txt > temp
cat header>france.asc
cat temp | awk '{printf "%10.6f%10.6f%10.6f%10.6f\n", $1,-$2,1,1}' >>france.asc
echo 'END      3.33e+032'>>france.asc

Here is the ASCII file thus created (gzipped) and I converted it into binary NTv2 format (as used by PROJ.4) using the GDAit Windows software:

france.gsb

To use the above file in GRASS, copy it to the $(GISBASE)/etc/nad directory, and add the line 'nadgrids: france.gsb' to your PROJ_INFO file in the PERMANENT directory of your ntf/clark80 location. Make sure there are no other towgs84, nadgrids or dx, dy or dz lines in the PROJ_INFO file.

Some final tests:
(output from the two commands should be the same)

echo '-2.9 46' | \
cs2cs -f %.3f +proj=longlat +ellps=GRS80 +to +proj=geocent +ellps=GRS80 | \
awk '{printf "%.3f\t%.3f\t%.3f\n",$1-(-168.882),$2-(-62.258),$3-(320.741)}' | \
cs2cs +proj=geocent +ellps=clrk80 +to +proj=longlat +ellps=clrk80


echo '-2.9 46' | \
cs2cs +proj=longlat +datum=WGS84 +to +proj=longlat +ellps=clrk80 \
+nadgrids=./france.gsb

Also this page in the PROJ.4 mailing list archive gives some hints on setting up the NTF Lambert projections used in France with PROJ.4. Just replace '+towgs84=-168,-60,+320' with '+nadgrids=./france.gsb' to make use of the higher accuracy datum shift file.

You may also compare the results with those from the Circé 2000 co-ordinate convertor available from the IGN Téléchargement page mentioned above. When we get down to 8 or 9 decimal places the results don't agree; I think this is because of the different interpolation methods used. If you want your co-ordinate conversions to be the official very correct values, you would need to implement support for reading the French grid-shift file directly in PROJ.4 or whatever co-ordinate software you use. The file I have provided here is just to make it easy and convenient to get most of the accuracy without having to write any new software. The main motivation is of course use in GRASS, which has the capability to re-project raster and vector maps as well as lists of sites.

Shift

Accuracy