Amersfoort / RD New

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

Amersfoort / RD New

Roger Oberholtzer-2
I am trying to get the Amersfoort / RD New projection working with
proj4. The users (a road authority in the Netherlands) claim to use
this definition with ESRI:

PROJCS["RD_New",GEOGCS["GCS_Amersfoort",DATUM["D_Amersfoort",
SPHEROID["Bessel_1841",6377397.155,299.1528128]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Double_Stereographic"],
PARAMETER["False_Easting",155000.0],
PARAMETER["False_Northing",463000.0],
PARAMETER["Central_Meridian",5.38763888888889],
PARAMETER["Scale_Factor",0.9999079],
PARAMETER["Latitude_Of_Origin",52.15616055555555],UNIT["Meter",1.0]]

I have started with the definition at
http://spatialreference.org/ref/epsg/amersfoort-rd-new/proj4/. The
projections from that are 100 meters from where they should be. I have
seen claims that the Spatial Reference definition is incorrect, and
that a better one is (it adds the +towgs84):

“+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
+towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725
+units=m +no_defs”

Still no joy. My source values are WGS84 lat/long from a GPS receiver.

Anyone have the correct proj4 definition for this?


--
Roger Oberholtzer
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Sebastiaan Couwenberg
On 2016-11-01 12:33, Roger Oberholtzer wrote:

> I am trying to get the Amersfoort / RD New projection working with
> proj4. The users (a road authority in the Netherlands) claim to use
> this definition with ESRI:
>
> PROJCS["RD_New",GEOGCS["GCS_Amersfoort",DATUM["D_Amersfoort",
> SPHEROID["Bessel_1841",6377397.155,299.1528128]],
> PRIMEM["Greenwich",0.0],
> UNIT["Degree",0.0174532925199433]],
> PROJECTION["Double_Stereographic"],
> PARAMETER["False_Easting",155000.0],
> PARAMETER["False_Northing",463000.0],
> PARAMETER["Central_Meridian",5.38763888888889],
> PARAMETER["Scale_Factor",0.9999079],
> PARAMETER["Latitude_Of_Origin",52.15616055555555],UNIT["Meter",1.0]]
>
> I have started with the definition at
> http://spatialreference.org/ref/epsg/amersfoort-rd-new/proj4/. The
> projections from that are 100 meters from where they should be. I have
> seen claims that the Spatial Reference definition is incorrect, and
> that a better one is (it adds the +towgs84):
>
> “+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
> +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725
> +units=m +no_defs”
>
> Still no joy. My source values are WGS84 lat/long from a GPS receiver.
>
> Anyone have the correct proj4 definition for this?

The +towgs84 parameters as published by the Dutch Kadaster are:

  +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.877402,4.0725

These are used in the epsg file from PROJ.4 4.9.3, see:

  https://github.com/osgeonl/rdprojectie/blob/master/doc/ntv2grid.md

Alternatively you can use the RDNAP grid shift files also published by
Kadaster. If you're on Debian/Ubuntu you can install the proj-rdnap
package (available in jessie-backports/non-free and Ubuntu wily and
later), otherwise you can find the download link in the PROJ.4
documentation at:

  http://proj4.org/grids.html#netherlands
  http://www.kadaster.nl/transformatie-van-coordinaten

Use the link in the "Benadering transformatieprocedure" section to go to
the request package and fill in the form. You'll be sent a download link
by email.

We've recently had a mini-seminar about the correct usage of RD in open
source software, and started collecting reference information resulting
from that in a git repository:

  https://github.com/osgeonl/rdprojectie

You'll find more information in the linked posts on the Dutch
mailinglist.

Kind Regards,

Bas
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Jelmer Baas
As an adition:

If you're trying to match the MapInfo (12.5 and earlier)  projection system, use this:
+proj=stere +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs

This is something I've been struggling with for a long time, and with help from a dev from within MapInfo, finally managed to find out.
Would appreciate it if someone could make a note of this, somewhere :)

Regards,
Jelmer Baas
Speer IT B.V.


-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of Bas Couwenberg
Sent: dinsdag 1 november 2016 12:52
To: [hidden email]
Subject: Re: [Proj] Amersfoort / RD New

On 2016-11-01 12:33, Roger Oberholtzer wrote:

> I am trying to get the Amersfoort / RD New projection working with
> proj4. The users (a road authority in the Netherlands) claim to use
> this definition with ESRI:
>
> PROJCS["RD_New",GEOGCS["GCS_Amersfoort",DATUM["D_Amersfoort",
> SPHEROID["Bessel_1841",6377397.155,299.1528128]],
> PRIMEM["Greenwich",0.0],
> UNIT["Degree",0.0174532925199433]],
> PROJECTION["Double_Stereographic"],
> PARAMETER["False_Easting",155000.0],
> PARAMETER["False_Northing",463000.0],
> PARAMETER["Central_Meridian",5.38763888888889],
> PARAMETER["Scale_Factor",0.9999079],
> PARAMETER["Latitude_Of_Origin",52.15616055555555],UNIT["Meter",1.0]]
>
> I have started with the definition at
> http://spatialreference.org/ref/epsg/amersfoort-rd-new/proj4/. The
> projections from that are 100 meters from where they should be. I have
> seen claims that the Spatial Reference definition is incorrect, and
> that a better one is (it adds the +towgs84):
>
> “+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
> +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725
> +units=m +no_defs”
>
> Still no joy. My source values are WGS84 lat/long from a GPS receiver.
>
> Anyone have the correct proj4 definition for this?

The +towgs84 parameters as published by the Dutch Kadaster are:

  +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.877402,4.0725

These are used in the epsg file from PROJ.4 4.9.3, see:

  https://github.com/osgeonl/rdprojectie/blob/master/doc/ntv2grid.md

Alternatively you can use the RDNAP grid shift files also published by
Kadaster. If you're on Debian/Ubuntu you can install the proj-rdnap
package (available in jessie-backports/non-free and Ubuntu wily and
later), otherwise you can find the download link in the PROJ.4
documentation at:

  http://proj4.org/grids.html#netherlands
  http://www.kadaster.nl/transformatie-van-coordinaten

Use the link in the "Benadering transformatieprocedure" section to go to
the request package and fill in the form. You'll be sent a download link
by email.

We've recently had a mini-seminar about the correct usage of RD in open
source software, and started collecting reference information resulting
from that in a git repository:

  https://github.com/osgeonl/rdprojectie

You'll find more information in the linked posts on the Dutch
mailinglist.

Kind Regards,

Bas
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Roger Oberholtzer-2
In reply to this post by Sebastiaan Couwenberg
I downloaded the grid shift file from http://www.kadaster.nl/transformatie-van-coordinaten

I am not using it yet. I took one of the sample points from the Use doc that came with it to see what I got compared to the expected. For example:

Zuid-Limburg:
        Ref:    Latitude:   50.7925849160 Longitude:    5.7737955480 Altitude:    245.948
        Ref:    Easting:       182260.450 Northing:       311480.670 Altitude:    200.000
        Proj4:  Easting:       182225.909 Northing:       311388.163 Altitude:    245.948    98.746 meters from Ref

The reference values are ETRS89. But I cannot think that accounts to 100 meters of difference. My proj spec is:

     +proj=sterea
     +lat_0=52.15616055555555
     +lon_0=5.38763888888889
     +k=0.9999079
     +x_0=155000
     +y_0=463000
     +ellps=bessel
     +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725
     +units=m
     +no_defs

Doing this via the C API as I do for many other projections. These are passed to pj_init_plus.

Would the grid shift make such a large correction?

On Tue, Nov 1, 2016 at 12:51 PM, Bas Couwenberg <[hidden email]> wrote:
On 2016-11-01 12:33, Roger Oberholtzer wrote:
> I am trying to get the Amersfoort / RD New projection working with
> proj4. The users (a road authority in the Netherlands) claim to use
> this definition with ESRI:
>
> PROJCS["RD_New",GEOGCS["GCS_Amersfoort",DATUM["D_Amersfoort",
> SPHEROID["Bessel_1841",6377397.155,299.1528128]],
> PRIMEM["Greenwich",0.0],
> UNIT["Degree",0.0174532925199433]],
> PROJECTION["Double_Stereographic"],
> PARAMETER["False_Easting",155000.0],
> PARAMETER["False_Northing",463000.0],
> PARAMETER["Central_Meridian",5.38763888888889],
> PARAMETER["Scale_Factor",0.9999079],
> PARAMETER["Latitude_Of_Origin",52.15616055555555],UNIT["Meter",1.0]]
>
> I have started with the definition at
> http://spatialreference.org/ref/epsg/amersfoort-rd-new/proj4/. The
> projections from that are 100 meters from where they should be. I have
> seen claims that the Spatial Reference definition is incorrect, and
> that a better one is (it adds the +towgs84):
>
> “+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
> +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725
> +units=m +no_defs”
>
> Still no joy. My source values are WGS84 lat/long from a GPS receiver.
>
> Anyone have the correct proj4 definition for this?

The +towgs84 parameters as published by the Dutch Kadaster are:

  +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.877402,4.0725

These are used in the epsg file from PROJ.4 4.9.3, see:

  https://github.com/osgeonl/rdprojectie/blob/master/doc/ntv2grid.md

Alternatively you can use the RDNAP grid shift files also published by
Kadaster. If you're on Debian/Ubuntu you can install the proj-rdnap
package (available in jessie-backports/non-free and Ubuntu wily and
later), otherwise you can find the download link in the PROJ.4
documentation at:

  http://proj4.org/grids.html#netherlands
  http://www.kadaster.nl/transformatie-van-coordinaten

Use the link in the "Benadering transformatieprocedure" section to go to
the request package and fill in the form. You'll be sent a download link
by email.

We've recently had a mini-seminar about the correct usage of RD in open
source software, and started collecting reference information resulting
from that in a git repository:

  https://github.com/osgeonl/rdprojectie

You'll find more information in the linked posts on the Dutch
mailinglist.

Kind Regards,

Bas
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj



--
Roger Oberholtzer

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Sebastiaan Couwenberg
On 2016-11-01 15:05, Roger Oberholtzer wrote:

> I downloaded the grid shift file from
> http://www.kadaster.nl/transformatie-van-coordinaten
>
> I am not using it yet. I took one of the sample points from the Use doc
> that came with it to see what I got compared to the expected. For
> example:
>
> Zuid-Limburg:
> Ref:   Latitude: 50.7925849160 Longitude: 5.7737955480 Altitude:
> 245.948
> Ref:   Easting:     182260.450 Northing:    311480.670 Altitude:
> 200.000
> Proj4: Easting:     182225.909 Northing:    311388.163 Altitude:
> 245.948
> 98.746 meters from Ref
>
> The reference values are ETRS89. But I cannot think that accounts to
> 100
> meters of difference. My proj spec is:
>
>      +proj=sterea
>      +lat_0=52.15616055555555
>      +lon_0=5.38763888888889
>      +k=0.9999079
>      +x_0=155000
>      +y_0=463000
>      +ellps=bessel
>
> +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725
>      +units=m
>      +no_defs
>
> Doing this via the C API as I do for many other projections. These are
> passed to pj_init_plus.
>
> Would the grid shift make such a large correction?

 From the proj-rdnap test script (using PROJ.4 4.9.3):

  Test:   05 Zuid-Limburg
  Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
  Input:  50.792584908 5.773795547 174.9478
  Output: 182260.4500     311480.6701 129.0000
  Expect: 182260.4500     311480.6700 129.000
  Test OK: From ETRS89 to RD/NAP - 05 Zuid-Limburg (Not identical, but
within margin)

The above uses the test values from "Use of RDTRANS2008 and
NAPTRANS2008.pdf" included in the NTv2.zip, your values are different,
not sure where you got those.

The rdnap configuration for the above test contains:

  # RDNAP with NTv2 and VDatum
  <rdnap> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
+nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m +no_defs
<>

Try using the grid shift files for your test, it should be more accurate
(at ground level) than using +towgs84 values.

Kind Regards,

Bas
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Roger Oberholtzer-2
I have updated to 4.9.3. My rdnap file is identical to yours:

# RDNAP with NTv2 and VDatum
<rdnap> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
+nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m +no_defs
<>

If I run this:
cs2cs  -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f

And provide these locations:

50.792584908 5.773795547 174.9478

I get this:

182225.9086     311388.1617 128.9976

Not the values you listed:

182260.4500     311480.6701 129.0000


I'm confused.



On Tue, Nov 1, 2016 at 3:28 PM, Bas Couwenberg <[hidden email]> wrote:
On 2016-11-01 15:05, Roger Oberholtzer wrote:
> I downloaded the grid shift file from
> http://www.kadaster.nl/transformatie-van-coordinaten
>
> I am not using it yet. I took one of the sample points from the Use doc
> that came with it to see what I got compared to the expected. For
> example:
>
> Zuid-Limburg:
> Ref:   Latitude: 50.7925849160 Longitude: 5.<a href="tel:7737955480" value="+17737955480">7737955480 Altitude:
> 245.948
> Ref:   Easting:     182260.450 Northing:    311480.670 Altitude:
> 200.000
> Proj4: Easting:     182225.909 Northing:    311388.163 Altitude:
> 245.948
> 98.746 meters from Ref
>
> The reference values are ETRS89. But I cannot think that accounts to
> 100
> meters of difference. My proj spec is:
>
>      +proj=sterea
>      +lat_0=52.15616055555555
>      +lon_0=5.38763888888889
>      +k=0.9999079
>      +x_0=155000
>      +y_0=463000
>      +ellps=bessel
>
> +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725
>      +units=m
>      +no_defs
>
> Doing this via the C API as I do for many other projections. These are
> passed to pj_init_plus.
>
> Would the grid shift make such a large correction?

 From the proj-rdnap test script (using PROJ.4 4.9.3):

  Test:   05 Zuid-Limburg
  Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
  Input:  50.792584908 5.773795547 174.9478
  Output: 182260.4500     311480.6701 129.0000
  Expect: 182260.4500     311480.6700 129.000
  Test OK: From ETRS89 to RD/NAP - 05 Zuid-Limburg (Not identical, but
within margin)

The above uses the test values from "Use of RDTRANS2008 and
NAPTRANS2008.pdf" included in the NTv2.zip, your values are different,
not sure where you got those.

The rdnap configuration for the above test contains:

  # RDNAP with NTv2 and VDatum
  <rdnap> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
+nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m +no_defs
<>

Try using the grid shift files for your test, it should be more accurate
(at ground level) than using +towgs84 values.

Kind Regards,

Bas
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj



--
Roger Oberholtzer

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Sebastiaan Couwenberg
On 2016-11-02 09:36, Roger Oberholtzer wrote:

> I have updated to 4.9.3. My rdnap file is identical to yours:
>
> # RDNAP with NTv2 and VDatum
> <rdnap> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
> +nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m
> +no_defs
> <>
>
> If I run this:
> cs2cs  -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
>
> And provide these locations:
>
> 50.792584908 5.773795547 174.9478
>
> I get this:
>
> 182225.9086     311388.1617 128.9976
>
> Not the values you listed:
>
> 182260.4500     311480.6701 129.0000
>
>
> I'm confused.

Does your epsg file use +towgs84 values for EPSG:4258?

  # ETRS89
  <4258> +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  <>

I guess you're on Windows, which uses a different libc that may also
explain the divergence.

> On Tue, Nov 1, 2016 at 3:28 PM, Bas Couwenberg <[hidden email]>
> wrote:
>
>> On 2016-11-01 15:05, Roger Oberholtzer wrote:
>> > I downloaded the grid shift file from
>> > http://www.kadaster.nl/transformatie-van-coordinaten
>> >
>> > I am not using it yet. I took one of the sample points from the Use doc
>> > that came with it to see what I got compared to the expected. For
>> > example:
>> >
>> > Zuid-Limburg:
>> > Ref:   Latitude: 50.7925849160 Longitude: 5.7737955480 Altitude:
>> > 245.948
>> > Ref:   Easting:     182260.450 Northing:    311480.670 Altitude:
>> > 200.000
>> > Proj4: Easting:     182225.909 Northing:    311388.163 Altitude:
>> > 245.948
>> > 98.746 meters from Ref
>> >
>> > The reference values are ETRS89. But I cannot think that accounts to
>> > 100
>> > meters of difference. My proj spec is:
>> >
>> >      +proj=sterea
>> >      +lat_0=52.15616055555555
>> >      +lon_0=5.38763888888889
>> >      +k=0.9999079
>> >      +x_0=155000
>> >      +y_0=463000
>> >      +ellps=bessel
>> >
>> > +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.
>> 343987817378283,-1.87740163998045,4.0725
>> >      +units=m
>> >      +no_defs
>> >
>> > Doing this via the C API as I do for many other projections. These are
>> > passed to pj_init_plus.
>> >
>> > Would the grid shift make such a large correction?
>>
>>  From the proj-rdnap test script (using PROJ.4 4.9.3):
>>
>>   Test:   05 Zuid-Limburg
>>   Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
>>   Input:  50.792584908 5.773795547 174.9478
>>   Output: 182260.4500     311480.6701 129.0000
>>   Expect: 182260.4500     311480.6700 129.000
>>   Test OK: From ETRS89 to RD/NAP - 05 Zuid-Limburg (Not identical, but
>> within margin)
>>
>> The above uses the test values from "Use of RDTRANS2008 and
>> NAPTRANS2008.pdf" included in the NTv2.zip, your values are different,
>> not sure where you got those.
>>
>> The rdnap configuration for the above test contains:
>>
>>   # RDNAP with NTv2 and VDatum
>>   <rdnap> +proj=sterea +lat_0=52.15616055555555
>> +lon_0=5.38763888888889
>> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
>> +nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m
>> +no_defs
>> <>
>>
>> Try using the grid shift files for your test, it should be more
>> accurate
>> (at ground level) than using +towgs84 values.
>>
>> Kind Regards,
>>
>> Bas
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Roger Oberholtzer-2


On Wed, Nov 2, 2016 at 9:48 AM, Bas Couwenberg <[hidden email]> wrote:
On 2016-11-02 09:36, Roger Oberholtzer wrote:
> I have updated to 4.9.3. My rdnap file is identical to yours:
>
> # RDNAP with NTv2 and VDatum
> <rdnap> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
> +nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m
> +no_defs
> <>
>
> If I run this:
> cs2cs  -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
>
> And provide these locations:
>
> 50.792584908 5.773795547 174.9478
>
> I get this:
>
> 182225.9086     311388.1617 128.9976
>
> Not the values you listed:
>
> 182260.4500     311480.6701 129.0000
>
>
> I'm confused.

Does your epsg file use +towgs84 values for EPSG:4258?

  # ETRS89
  <4258> +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  <>

I guess you're on Windows, which uses a different libc that may also
explain the divergence.

I'm on Linux. I see that when I updated my libproj, a new epsg file was installed. I don't usually use that file. I did here because I wanted to duplicate your exact command. 4258 is defined as:
 
<4258> +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  <>

I'm curious why you use both the definition in the epsg file, as well as the one in rdnap. Why not everything in the rdnap file? What is your definition for 4258?

If I change my command to:

cs2cs  -r  +to +init=rdnap:rdnap

I see that the grid shift file is now read. But my computed values are still:

182225.9086     311388.1617 128.9976




--
Roger Oberholtzer

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Sebastiaan Couwenberg
On 11/02/2016 11:06 AM, Roger Oberholtzer wrote:

> On Wed, Nov 2, 2016 at 9:48 AM, Bas Couwenberg <[hidden email]> wrote:
>
>> On 2016-11-02 09:36, Roger Oberholtzer wrote:
>>> I have updated to 4.9.3. My rdnap file is identical to yours:
>>>
>>> # RDNAP with NTv2 and VDatum
>>> <rdnap> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
>>> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
>>> +nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m
>>> +no_defs
>>> <>
>>>
>>> If I run this:
>>> cs2cs  -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
>>>
>>> And provide these locations:
>>>
>>> 50.792584908 5.773795547 174.9478
>>>
>>> I get this:
>>>
>>> 182225.9086     311388.1617 128.9976
>>>
>>> Not the values you listed:
>>>
>>> 182260.4500     311480.6701 129.0000
>>>
>>>
>>> I'm confused.
>>
>> Does your epsg file use +towgs84 values for EPSG:4258?
>>
>>   # ETRS89
>>   <4258> +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  <>
>>
>> I guess you're on Windows, which uses a different libc that may also
>> explain the divergence.
>>
>
> I'm on Linux. I see that when I updated my libproj, a new epsg file was
> installed. I don't usually use that file. I did here because I wanted to
> duplicate your exact command. 4258 is defined as:
>
> <4258> +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  <>
>
> I'm curious why you use both the definition in the epsg file, as well as
> the one in rdnap. Why not everything in the rdnap file? What is your
> definition for 4258?

My definition for EPSG:4258 is the one from PROJ.4 4.9.3, the one below
"# ETRS89" in the quote above.

Since the EPSG:4258 definition is already available in the epsg file,
there is no need to duplicate it in the rdnap file.

> If I change my command to:
>
> cs2cs  -r  +to +init=rdnap:rdnap
>
> I see that the grid shift file is now read. But my computed values are
> still:
>
> 182225.9086     311388.1617 128.9976

We'll need more information about your environment to find out the cause
of the difference. I'm starting to suspect that the old PROJ.4
installation is still being used, not the new version you installed
elsewhere.

I cannot reproduce your problem, both PROJ.4 4.8.0 on Debian stable &
4.9.3 on Debian unstable produce the same result for me:

 $ lsb_release -ds
 Debian GNU/Linux 8.6 (jessie)
 $ apt-cache show proj-bin | grep Version
 Version: 4.8.0-5
 $ echo "50.792584908 5.773795547 174.9478" \
 | PROJ_DEBUG=3 cs2cs  -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
 pj_open_lib(epsg): call fopen(/usr/share/proj/epsg) - succeeded

 pj_open_lib(rdnap): call fopen(/usr/share/proj/rdnap) - succeeded

 pj_open_lib(rdtrans2008.gsb): call
fopen(/usr/share/proj/rdtrans2008.gsb) - succeeded

 NTv2 NL_ALL   63x65: LL=(2.5,50.5) UR=(7.66666667,55.8333333)

 NTv2 NL_LAND  561x421: LL=(3,50.5) UR=(7.66666667,54)

 NTv2 - loading grid NL_LAND
 pj_open_lib(rdtrans2008.gsb): call
fopen(/usr/share/proj/rdtrans2008.gsb) - succeeded

 pj_apply_gridshift(): used NL_LAND
 pj_open_lib(naptrans2008.gtx): call
fopen(/usr/share/proj/naptrans2008.gtx) - succeeded

 GTX 311x641: LL=(2.5,50.5) UR=(7.66666667,55.8333333)
 pj_open_lib(naptrans2008.gtx): call
fopen(/usr/share/proj/naptrans2008.gtx) - succeeded

 pj_apply_gridshift(): used GTX Vertical Grid Shift File
 182260.4500     311480.6701 129.0000


 $ lsb_release -ds
 Debian GNU/Linux unstable (sid)
 $ apt-cache show proj-bin | grep Version
 Version: 4.9.3-1
 $ echo "50.792584908 5.773795547 174.9478" \
 | PROJ_DEBUG=3 cs2cs  -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
 pj_open_lib(epsg): call fopen(/usr/share/proj/epsg) - succeeded

 pj_open_lib(rdnap): call fopen(/usr/share/proj/rdnap) - succeeded

 pj_open_lib(rdtrans2008.gsb): call
fopen(/usr/share/proj/rdtrans2008.gsb) - succeeded

 NTv2 NL_ALL   63x65: LL=(2.5,50.5) UR=(7.66666667,55.8333333)

 NTv2 NL_LAND  561x421: LL=(3,50.5) UR=(7.66666667,54)

 NTv2 - loading grid NL_LAND
 pj_open_lib(rdtrans2008.gsb): call
fopen(/usr/share/proj/rdtrans2008.gsb) - succeeded

 pj_apply_gridshift(): used NL_LAND
 pj_open_lib(naptrans2008.gtx): call
fopen(/usr/share/proj/naptrans2008.gtx) - succeeded

 GTX 311x641: LL=(2.5,50.5) UR=(7.66666667,55.8333333)
 pj_open_lib(naptrans2008.gtx): call
fopen(/usr/share/proj/naptrans2008.gtx) - succeeded

 pj_apply_gridshift(): used GTX Vertical Grid Shift File
 182260.4500     311480.6701 129.0000


Kind Regards,

Bas

--
 GPG Key ID: 4096R/6750F10AE88D4AF1
Fingerprint: 8182 DE41 7056 408D 6146  50D1 6750 F10A E88D 4AF1
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Roger Oberholtzer-2
It could be that the old epsg file was still being used. The old (non-working) definition of 4258 was:

     <4258> +proj=longlat +ellps=GRS80 +no_defs  <>

and the new (working) is:

     <4258> +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  <>

I think the problem was the missing +towgs84=0,0,0,0,0,0,0. Without that, the grid shift file was not used. And probably other things as well.

I want to continue this without reference to files as they are too variable to be reliable. Instead, I have moved everything to the command line so there is no confusion about what is in each system's files.

The command:

     echo "50.792584908 5.773795547 174.9478" | \
         cs2cs -r +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs \
              +to +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 \
                +x_0=155000 +y_0=463000 +ellps=bessel +nadgrids=rdtrans2008.gsb \
                +geoidgrids=naptrans2008.gtx +units=m +no_defs -f %.4f

Now correctly provides:

182260.4500     311480.6701 129.0000

I have 'translated' this to the C API and all is working great. I have been testing with the reference points that came with the grid file.

Texel:
        Ref:    Latitude:   53.1607530420 Longitude:    4.8247619120 Altitude:     42.861
        Ref:    Easting:       117380.120 Northing:       575040.340 Altitude:      1.000
        Calc:   Easting:       117380.120 Northing:       575040.340 Altitude:      1.000    0.000 meters from Ref
        Inv:    Latitude:   53.1607530431 Longitude:    4.8247619081 Altitude:   2455.773

Noord-Groningen:
        Ref:    Latitude:   53.4194820500 Longitude:    6.7767266740 Altitude:     42.359
        Ref:    Easting:       247380.560 Northing:       604580.780 Altitude:      2.000
        Calc:   Easting:       247380.560 Northing:       604580.780 Altitude:      2.000    0.000 meters from Ref
        Inv:    Latitude:   53.4194820513 Longitude:    6.7767266740 Altitude:   2426.971

Amersfoort:
        Ref:    Latitude:   52.1551728970 Longitude:    5.3872036570 Altitude:     43.255
        Ref:    Easting:       155000.000 Northing:       463000.000 Altitude:      0.000
        Calc:   Easting:       155000.000 Northing:       463000.000 Altitude:      0.000    0.000 meters from Ref
        Inv:    Latitude:   52.1551728982 Longitude:    5.3872036577 Altitude:   2478.331

Amersfoort 100m:
        Ref:    Latitude:   52.1551729100 Longitude:    5.3872036580 Altitude:    143.255
        Ref:    Easting:       155000.000 Northing:       463000.000 Altitude:    100.000
        Calc:   Easting:       155000.000 Northing:       463000.001 Altitude:    100.000    0.001 meters from Ref
        Inv:    Latitude:   52.1551728982 Longitude:    5.3872036577 Altitude:   8207.909

Zeeuws-Vlaanderen:
        Ref:    Latitude:   51.3686071520 Longitude:    3.3975885950 Altitude:     47.402
        Ref:    Easting:        16460.910 Northing:       377380.230 Altitude:      3.000
        Calc:   Easting:        16460.910 Northing:       377380.230 Altitude:      3.000    0.000 meters from Ref
        Inv:    Latitude:   51.3686071526 Longitude:    3.3975885974 Altitude:   2715.954

Zuid-Limburg:
        Ref:    Latitude:   50.7925849160 Longitude:    5.7737955480 Altitude:    245.948
        Ref:    Easting:       182260.450 Northing:       311480.670 Altitude:    200.000
        Calc:   Easting:       182260.450 Northing:       311480.671 Altitude:    200.000    0.001 meters from Ref
        Inv:    Latitude:   50.7925849065 Longitude:    5.7737955467 Altitude:  14091.766

Maasvlakte:
        Ref:    Latitude:   51.9473938980 Longitude:    4.0728871010 Altitude:     47.597
        Ref:    Easting:        64640.890 Northing:       440700.010 Altitude:      4.000
        Calc:   Easting:        64640.890 Northing:       440700.010 Altitude:      4.000    0.000 meters from Ref

        Inv:    Latitude:   51.9473938985 Longitude:    4.0728871025 Altitude:   2727.102

The Ref values are the sample points. The Calc is what I now get from proj (yeah!). The Inv is what I get when I do an inverse projection of my projected point. It is rather close to the original Ref. Except the Altitude. I am using the shift file for that as well. Maybe it does not allow an inverse on the Altitude? It is probably not a big issue for me. More of a curiosity.

I am also curious about some additional sample points that came with the grid shift data.

outside:
        Ref:    Latitude:   48.8430302100 Longitude:    8.7232602350 Altitude:     52.029
        Ref:    Easting:       400000.230 Northing:       100000.450 Altitude:      5.000

no_rd&geoid:
        Ref:    Latitude:   50.6874203920 Longitude:    4.6089718130 Altitude:     51.611
        Ref:    Easting:       100000.670 Northing:       300000.890 Altitude:      6.000

no_geoid:
        Ref:    Latitude:   51.1368251970 Longitude:    4.6013753610 Altitude:     50.967
        Ref:    Easting:       100000.670 Northing:       350000.890 Altitude:      6.000

no_rd:
        Ref:    Latitude:   52.4824408390 Longitude:    4.2684038890 Altitude:     49.944
        Ref:    Easting:        79000.010 Northing:       500000.230 Altitude:      7.000

edge_rd:
        Ref:    Latitude:   51.0039765320 Longitude:    3.8912478300 Altitude:     52.743
        Ref:    Easting:        50000.450 Northing:       335999.670 Altitude:      8.000


All of these fail. Since I expect our data to be within the area of the grid shift file, I guess I do not need to worry about these. Out of curiosity, how do others deal with such points?

--
Roger Oberholtzer

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Amersfoort / RD New

Sebastiaan Couwenberg
On 2016-11-03 09:15, Roger Oberholtzer wrote:

> It could be that the old epsg file was still being used. The old
> (non-working) definition of 4258 was:
>
>      <4258> +proj=longlat +ellps=GRS80 +no_defs  <>
>
> and the new (working) is:
>
>      <4258> +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  
> <>
>
> I think the problem was the missing +towgs84=0,0,0,0,0,0,0. Without
> that,
> the grid shift file was not used. And probably other things as well.

Yes, the missing +towgs84 values are the cause, I experienced this too
with the epsg file from PROJ.4 4.7.0 which doesn't have the +towgs84
values.

See: https://lists.osgeo.org/pipermail/dutch/2014-November/001014.html

> The Ref values are the sample points. The Calc is what I now get from
> proj
> (yeah!). The Inv is what I get when I do an inverse projection of my
> projected point. It is rather close to the original Ref. Except the
> Altitude. I am using the shift file for that as well. Maybe it does not
> allow an inverse on the Altitude? It is probably not a big issue for
> me.
> More of a curiosity.

The documentation included along with the grid shift files documents the
altitude limitation:

"
  1) The rdtrans2008 NTv2-grid can only give identical results to
RDNAPTRANS TM 2008 within 1
     millimeter at ground level onshore and at mean seal level offshore.
The horizontal deviation is
     approximately 1 millimeter per 50 meter height difference from
ground level or mean sea level.
  2) An exception to 1) is the border of the RDNAPTRANS TM 2008
correction grid. Transformation
     results within cells of the rdtrans2008 NTv2-grid that are
intersected by the border of the
     RDNAPTRANS TM 2008 correction grid can result in deviations of up to
20 centimeter.
  3) The naptrans2008 VDatum-grid cannot be used to determine deflections
of the vertical. For
     this the NLGEO2004 geoid model has to be used.
  4) The naptrans2008 VDatum-grid is referenced to the Bessel-1841
ellipsoid and cannot be used
     stand-alone, it has to be used in combination with the rdtrans2008
NTv2-grid.
"

https://anonscm.debian.org/cgit/pkg-grass/proj-rdnap.git/plain/Use%20of%20RDTRANS2008%20and%20NAPTRANS2008.pdf

> I am also curious about some additional sample points that came with
> the
> grid shift data.
>
> [...]
>
> All of these fail. Since I expect our data to be within the area of the
> grid shift file, I guess I do not need to worry about these. Out of
> curiosity, how do others deal with such points?

This also noted in the documentation for the grid shift files:

"
   Points 07 - 10 are outside the region where interpolation between
   either the NLGEO2004 geoid or the RD correction grid points is
   possible. RD is defined only within the region enclosed by the
   following points (in RD), outside this region RD coordinates can
   be computed, but the output should be handled with care.

   Corners of the validity region for RD:
   ┌────────┬────────┐
   │  x (m) │  y (m) │
   ├────────┼────────┤
   │ 140000 │ 630000 │
   │ 100000 │ 600000 │
   │  80000 │ 500000 │
   │  -8000 │ 390000 │
   │  -8000 │ 335000 │
   │ 100000 │ 335000 │
   │ 160000 │ 288000 │
   │ 220000 │ 288000 │
   │ 301000 │ 450000 │
   │ 301000 │ 615000 │
   │ 260000 │ 630000 │
   └────────┴────────┘
"

With the huge values change in PROJ.4 4.9.2 some of these values could
be calculated correctly, but caused most other tests to fail, so that
was reverted. See the discussion I linked before:

  http://lists.maptools.org/pipermail/proj/2016-February/007327.html

In my test script I added an override to accept the '*' output from
cs2cs because these tests are expected to fail.

Kind Regards,

Bas
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj