[gdal-dev] How to set GDAL/cs2cs/proj to use NTv2 grid shift file when using EPSG codes?

Tue May 28 16:22:05 PDT 2013

I tested my modification to support the transformation from AGD66 to GDA94 in Australia using a NTv2 file. CRS details are below:

AGD66 => EPSG:4202 => +proj=longlat +ellps=aust_SA +towgs84=-117.808,-51.536,137.784,0.303,0.446,0.234,-0.29 +no_defs
                                       => +proj=longlat +ellps=aust_SA +nadgrids=A66_National_13.09.01.gsb +no_defs

GDA94 => EPSG:4283 => +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

To support a direct use of the NTv2 (.gsb) file when initializing with an EPSG code, 4202 in the epsg file is now:
	# AGD66
	<4202> +proj=longlat +datum=AGD66 +no_defs  <> 
and the datum is identified in pj_datums.c as:
	"AGD66", "nadgrids=A66_National_13.09.01.gsb", "aust_SA", "Australian Geodetic Datum 1966",

For a set of 4 points from the ICSM documentation, I tested this code modification of pyproj against osgeo/ogr, and cs2cs using two different proj strings to define the source CRS (one with nadgrids the other using towgs84) and the results are the following:
(Lat positive east, Lon positive north, Diff_ columns are the difference between the ICSM result and the calculated results in arc-seconds)

ICSM_inX	ICSM_outX	osgeoX		cs2csnadgridsX	cs2cstowgs84X	pyprojX	Diff_osgeoX	Diff_cs2csnadgridsX	Diff_cs2cstowgs84X	Diff_pyprojX
147.437368	147.438731	147.4387352	147.438731	147.438735	147.4387312	0.0152568	0.0			0.0144000		0.0008136
130.654978	130.656199	130.6562168	130.6562	130.656217	130.6561997	0.0639828	0.0036000		0.0648000		0.0024696
144.423552	144.424868	144.4248765	144.424868	144.424877	144.4248684	0.0307584	0.0			0.0324000		0.0016128
143.925176	143.926495	143.9265029	143.926496	143.926503	143.9264956	0.0284796	0.0036000		0.0288000		0.0020412

ICSM_inY	ICSM_outY	osgeoY		cs2csnadgridsY	cs2cstowgs84Y	pyprojY	Diff_osgeoY	Diff_cs2csnadgridsY	Diff_cs2cstowgs84Y	Diff_pyprojY
-42.806215	-42.804718	-42.80471228	-42.804718	-42.804712	-42.80471834	0.0205844	0.0			0.0216000		0.0012323
-18.027176	-18.025773	-18.02573847	-18.025773	-18.025738	-18.02577316	0.1243058	0.0			0.1260000		0.0005580
-37.952536	-37.951034	-37.95103726	-37.951034	-37.951037	-37.95103375	0.0117468	0.0			0.0108000		0.0009162
-37.654321	-37.652821	-37.65282602	-37.652821	-37.652826	-37.65282077	0.0180580	0.0			0.0180000		0.0008330

The smallest differences with the ICSM results are achieved with cs2cs using the nadgrids parameter and the modified pyproj.
I also used the osgeo python library and its output seems to be similar than cs2cs using the towgs84 parameter.

For my usage I do not need to get into sub-metric accuracy so I assume the cs2cs using the nadgrids parameter and the modified pyproj to be roughly the same.
I still wonder why there is a small difference as they should be using the same parameters. 
Any ideas? Due to float rounding being handled differently?

I am looking into expanding this to following datums as they are country specific:
 AGD 1966, AGD 1984, NAD83(HARN), Cape, Lisbon, Datum 73, DHDN, Corrego Alegre 1961, Corrego Alegre 1970-72, NZGD 1949, OSGB36.

I am now encountering some difficulties dealing with NTv2 files that are not covering the full extents of the area where a datum is being used.
Examples are:
- PENR2009.gsb for Spain converting from ED50 (4230) to ETRS89 (4258) but ED50 was used over other countries in Europe,
- SAD69_003.GSB and SAD96_003.GSB for Brazil converting from SAD (4291 / 4618) to SIRGAS2000 (4674) but SAD was used over other countries in South America.

I tried:
> cs2cs -v -f %.6f +proj=longlat +ellps=intl +nadgrids=@PENR2009.gsb +towgs84=-87,-98,-121,0,0,0,0 +no_defs 
>                           +to +proj=longlat +ellps=GRS80 
> +towgs84=0,0,0,0,0,0,0 +no_defs
But there is a mention that the towgs84 parameter was specified but not used and then the conversion outside of Spain is not done (not even with the towgs84 parameter)

How could we have proj reverting to the towgs84 parameter when a point falls outside of the area covered by the NTv2 file?

Thomas

-----Original Message-----
From: gdal-dev-bounces at lists.osgeo.org [mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of Thomas Campagne
Sent: May-22-13 11:43 AM
To: gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] How to set GDAL/cs2cs/proj to use NTv2 grid shift file when using EPSG codes?

Thanks André!

I will add a line for AGD66 datum and nadgrid file in pj_datums.c, and make the appropriate change in the epsg file , and compile from source. 

Is possible to default to the use of +towgs84 when there is no successful grids with +nadgrids?

For New-Zealand the line in pj_datums.c is:
"nzgd49",	"nadgrids=nzgd2kgrid0005.gsb", 	"intl",	 "New Zealand Geodetic Datum 1949"
And the epsg file contains:
# NZGD49
<4272> +proj=longlat +datum=nzgd49 +no_defs  <>

Could the following work if points are not located within the area covered by the NTv2 file?

pj_datums.c contains:
"nzgd49",	" nadgrids=@nzgd2kgrid0005.gsb ", 	"intl", "New Zealand Geodetic Datum 1949"
And the epsg file contains:
# NZGD49
<4272> +proj=longlat +datum=nzgd49 + towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993 +no_defs  <>

-----Original Message-----
From: gdal-dev-bounces at lists.osgeo.org [mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of Andre Joost
Sent: May-17-13 10:21 AM
To: gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] How to set GDAL/cs2cs/proj to use NTv2 grid shift file when using EPSG codes?

Am 17.05.2013 19:00, schrieb Thomas Campagne:

> When I run cs2cs -v +init=epsg:4267 I still do not understand how the
> +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat parameter is being
> generated as I cannot find it within the text/csv files.
>

EPSG:4267 expands to
+proj=longlat +datum=NAD27 +no_defs

and the nadgrids are in the "+datum=NAD27". This is hardcoded in <http://trac.osgeo.org/proj/browser/trunk/proj/src/pj_datums.c>
> "NAD27",    "nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat",
> 	                                                     "clrk66",
> 	                                "North_American_Datum_1927",

> Any idea about which files to edit to add the nadgrids parameter to 
> the proj string referenced by a given EPSG code (4202 for AGD66 as 
> example)?

You could add a line for AGD66 datum and nadgrid file in the file, and compile from source. And you have to change the share/proj/epsg file to use your datum istead of +towgs84.

Or try to insert the +nadgrids=<filename> for the AGD66 entry directly.

HTH,
André Joost

_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev