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 [email protected]
> +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: [email protected]
[mailto:[email protected]] On Behalf Of Thomas Campagne
Sent: May-22-13 11:43 AM
To: [email protected]
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", " [email protected] ", "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: [email protected]
[mailto:[email protected]] On Behalf Of Andre Joost
Sent: May-17-13 10:21 AM
To: [email protected]
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,@alaska,@ntv2_0.gsb,@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,@alaska,@ntv2_0.gsb,@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
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
