[Proj] proj and Dutch RD

Wed Oct 18 06:57:45 PDT 2006

On Wed, 2006-10-18 at 03:52 -0700, tvijlbrief wrote:
> No, it is not the same.
> 
>  "+proj=latlong +ellps=WGS84" will give incorrect results because there is
> no conversion
> from the Bessel to the WGS84 spheroid!

Isn't this part of the spec controlling what the source data is?
>From Bessel to wgs84? If my source data IS wgs84, why would I do this? 

I program from the C API where the specification of the source and the
specification of the destination are explicitly stated. Perhaps this is
why there is confusion on my part.

My code for RD is:

fromProj = pj_init_plus("+proj=latlong +ellps=WGS84");

toProj   = pj_init_plus("+proj=sterea +lat_0=52.15616055555555
				+lon_0=5.38763888888889
				+k=0.9999079 +x_0=155000 +y_0=463000
				+ellps=bessel +units=m +no_defs");

Then I use 

	lng = LONG * DEGREE_TO_RADIAN;
        lat = LAT  * DEGREE_TO_RADIAN;
        alt = 0; // We do not give this a height

        pj_transform(fromProj, toProj, 1, 0, &lng, &lat, &alt);

        local_easting = lng;
        local_northing = lat;

All with proper error checking.

> 
> You can compare the results by trying both:
> 
> tom at giant:~$ cs2cs +init=epsg:4326 +to +init=epsg:28992
> 6 53
> 196139.44       557179.10 -41.78
> tom at giant:~$ cs2cs +init=epsg:4326 +to +proj=sterea +lat_0=52.15616055555555
> +lon_0=5.38763888888889+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=WGS84
> +units=m +nodefs
> 6 53
> 196114.10       557077.21 0.00
> 
> The error is about 25 meters!!
> 
> Tom
> 
> 
> 
> Roger Oberholtzer wrote:
> > 
> > On Tue, 2006-10-17 at 02:35 -0700, tvijlbrief wrote:
> >> Much confusion, but the following definition will work! I submitted this
> >> as a patch for the epsg file to the libproj maintainers.
> >> 
> >> <28992> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 
> >> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel 
> >> +towgs84=565.040,49.910,465.840,-0.40939,0.35971,-1.86849,4.0772 +units=m 
> >> +no_defs  <>
> > 
> > Is the
> > "+towgs84=565.040,49.910,465.840,-0.40939,0.35971,-1.86849,4.0772" the
> > same as "+proj=latlong +ellps=WGS84" ?
> > 
> > I have an RD conversion and all mine is the same except for this.
> > 
> >> 
> >> For some reason Oscar van Vlijmen's conversion:
> >> 
> >> Transformation from Bessel 1841 to WGS84:
> >> dx=565.040; dy=49.910; dz=465.840; rx=0.40939; ry=-0.35971; rz=1.86849;
> >> ds=4.0772;
> >> 
> >> had the sign of rx,ry and rz reversed!
> >> 
> >> Tom
> >> 
> >> ====================
> >> 
> >> Roger Oberholtzer wrote:
> >> > 
> >> > On Thu, 2006-01-12 at 00:17 +0200, Martin Vermeer wrote:
> >> >> On Wed, Jan 11, 2006 at 12:43:13PM -0800, Eric Miller wrote:
> >> >> > I think Martin's comment is a little to brief.
> >> >> > 
> >> >> > Your source coordinate system is geodetic coordinates using the
> >> WGS84
> >> >> Datum (implies WGS84 ellipsoid)
> >> >> > In cs2cs terms:  +proj=latlong +datum=WGS84
> >> >> > 
> >> >> > Your destination coordinate system is RD (new)
> >> >> > In cs2cs terms: +proj=sterea +lat_0=52.15616055555555
> >> >> +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000
> >> >> +ellps=bessel +units=m +towgs84=593,26,478
> >> >> > 
> >> >> > I'm not sure about the signs for the towgs84 parameters (dx,dy,dz),
> >> >> maybe they're supposed to be negative.  Anyway, the +towgs84 tells
> >> cs2cs
> >> >> how to transform the WGS84 coordinates (in XYZ space) to whatever the
> >> >> output datum for RD is.  The +towgs84 or +datum for a different target
> >> >> coordinate systems will likely be different.  But, your source
> >> definition
> >> >> would stay the same.  If you have some reference coordinates, you can
> >> >> probably work out whether the towgs84 numbers are supposed to be
> >> positive
> >> >> or negative (provided those are good values).
> >> >> > 
> >> >> > Clear as mud?
> >> >> 
> >> >> This is precisely what I meant. Thanks!
> >> > 
> >> > But the 'fix' was to change the description of the source lat/longs
> >> from
> >> > 
> >> > 	+proj=latlong +ellps=WGS84
> >> > to
> >> > 	+proj=latlong +ellps=bessel
> >> > 
> >> > After the '+to' I always had +ellps=bessel to describe the desired
> >> > result. That was not changed.
> >> > 
> >> > Confusion here continues...
> >> > 
> >> >> 
> >> >> So play around with those +towgs84 parms, until you can reproduce your
> >> >> reference values. There must be decimetre precision shift values
> >> >> somewhere too, as I know the triangulation to be that precise. They
> >> must
> >> >> know at the cadastral survey in Apeldoorn. Here is their on-line
> >> >> transformer:
> >> >> 
> >> >> https://RDinfo.kadaster.nl/rd/transformator.html
> >> > 
> >> > I have looked at that. In summary, the conversion that does the best
> >> > (almost identical values) is:
> >> > 
> >> > +proj=latlong +ellps=bessel +to
> >> > +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
> >> > +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs
> >> > 
> >> > The one that was off by at least 60 meters was:
> >> > 
> >> > +proj=latlong +ellps=wgs84 +to
> >> > +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
> >> > +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs
> >> > 
> >> > -- 
> >> > Roger Oberholtzer
> >> > OPQ Systems AB
> >> > 
> >> > 
> >> > _______________________________________________
> >> > Proj mailing list
> >> > Proj at lists.maptools.org
> >> > http://lists.maptools.org/mailman/listinfo/proj
> >> > 
> >> > 
> >> 
> > -- 
> > Roger Oberholtzer
> > 
> > OPQ Systems AB
> > Ramböll Sverige AB
> > Kapellgränd 7
> > P.O. Box 4205
> > SE-102 65 Stockholm, Sweden
> > 
> > Tel: Int +46 8-615 60 20
> > Fax: Int +46 8-31 42 23
> > 
> > _______________________________________________
> > Proj mailing list
> > Proj at lists.maptools.org
> > http://lists.maptools.org/mailman/listinfo/proj
> > 
> 
-- 
Roger Oberholtzer

OPQ Systems AB
Ramböll Sverige AB
Kapellgränd 7
P.O. Box 4205
SE-102 65 Stockholm, Sweden

Tel: Int +46 8-615 60 20
Fax: Int +46 8-31 42 23