Problems with EOV for Hungary (1972 Egyseges Orszagos Vetuleti)

[Roger Oberholtzer]

I am trying to use proj (C API) to project WGS84 lat/longs in EOV for
Hungary (1972 Egyseges Orszagos Vetuleti). I am trying to match some
user data. The are using ArcGis with the following definition:

        PROJCS["HD_1972_Egyseges_Orszagos_Vetuleti",
        GEOGCS["GCS_Hungarian_1972",
        DATUM["D_Hungarian_1972",
        SPHEROID["GRS_1967",6378160,298.247167427]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.0174532925199432955]],
        PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
        PARAMETER["False_Easting",650000],
        PARAMETER["False_Northing",200000],
        PARAMETER["Scale_Factor",0.99993],
        PARAMETER["Azimuth",90],
        PARAMETER["Longitude_Of_Center",19.048571778],
        PARAMETER["Latitude_Of_Center",47.14439372222],
        UNIT["Meter",1]]

I have set this up as:

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

   to = pj_init_plus(
                "+proj=somerc "
                "+lat_0=47.14439372222222 "
                "+lon_0=19.04857177777778 "
                "+k=0.99993 "
                "+x_0=650000 "
                "+y_0=200000 "
                "+ellps=GRS67 "
                "+units=m "
                "+no_defs");

     pj_transform(to, from, 1, 0, &le, &ln, &la);

When I run some test data, I am always 80 some meters away. For example:

   Ref:    Latitude:   47.7656110000 Longitude:   18.9532900000
   Ref:    Easting:       642942.000 Northing:       269097.700
   Calc:   Easting:       642857.473 Northing:       269067.872

   89.636 meters from Ref

The two Ref. lines are the user data. The Calc. line is what I get with
the pj_init_plus.

I have tried other mercator proj= statements, but it seems to make no
difference. Does "Hotine_Oblique_Mercator_Azimuth_Center" have a
counterpart in proj?

Perhaps the GRS67 definition they use is not the same as proj's default?

The user claims that the WGS84 definition they use is:

        GEOGCS["GCS_WGS_1984",
        DATUM["D_WGS_1984",
        SPHEROID["WGS_1984",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]]

Isn't this the standard definition?

--
Roger Oberholtzer

OPQ Systems / Ramböll RST

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden

Office: Int [phone]
Mobile: Int [phone]

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

[Jean-Claude Repetto]

Roger Oberholtzer wrote:

>
>    to = pj_init_plus(
>                 "+proj=somerc "
>                 "+lat_0=47.14439372222222 "
>                 "+lon_0=19.04857177777778 "
>                 "+k=0.99993 "
>                 "+x_0=650000 "
>                 "+y_0=200000 "
>                 "+ellps=GRS67 "
>                 "+units=m "
>                 "+no_defs");
>
>      pj_transform(to, from, 1, 0, &le, &ln, &la);
>
> When I run some test data, I am always 80 some meters away. For example:

Hello,

You are missing the +towgs84 parameter, which is necessary if you want
PROJ.4 to perform a datum change.

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

[Mikael Rittri]

Jean-Claude wrote:
>
> Roger Oberholtzer wrote:
> > - - -
> > When I run some test data, I am always 80 some meters away.
>
>
> You are missing the +towgs84 parameter, which is necessary
> if you want PROJ.4 to perform a datum change.

Yes. I looked at www.epsg-registry.org, and searched on

   Type: Coordinate Transformation
   Area: Hungary

This gave a couple of 3-parameter and 7-parameter transforms
for HD72.  As usual, some of them were informally deprecated
in the comments field, but not formally DEPRECATED in the
technical sense.  Anyway, the best one seems to be

    EPSG:1448 "HD72 to WGS 84 (3)"

which has the same parameters as

    EPSG:1449 "HD72 to ETRS89 (2)",

namely
   
    X-axis translation   52.684   metre      
    Y-axis translation   -71.194  metre  
    Z-axis translation   -13.975  metre  
    X-axis rotation   0.312   arc-second  
    Y-axis rotation   0.1063  arc-second  
    Z-axis rotation   0.3729  arc-second  
    Scale difference  1.0191  parts per million  

However, there is a warning _above_ these parameters:

   "Care! Ensure application uses this method! See Guidance Note 7-1"

which means that you must scroll _down_ to the bottom line, where you find

   +Method [Coordinate Frame rotation]

and, since PROJ.4 uses the opposite sign convention for the rotation
angles (the Position Vector Transform), these angles must be negated.  
So, I think you should use

   +towgs84=52.684,-71.194,-13.975,-0.312,-0.1063,-0.3729,1.0191

to implement HD72.

Best regards,
--
Mikael Rittri
Carmenta AB
Göteborg
SWEDEN
www.carmenta.com

-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of Jean-Claude REPETTO
Sent: den 22 september 2009 11:57
To: PROJ.4 and general Projections Discussions; [hidden email]
Subject: Re: [Proj] Problems with EOV for Hungary (1972 Egyseges Orszagos Vetuleti)

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

[OvV_HN]

----- Original Message -----
From: "Roger Oberholtzer" <[hidden email]>
To: <[hidden email]>
Sent: Tuesday, September 22, 2009 11:02 AM
Subject: [Proj] Problems with EOV for Hungary (1972 Egyseges
OrszagosVetuleti)

>I am trying to use proj (C API) to project WGS84 lat/longs in EOV for
> Hungary (1972 Egyseges Orszagos Vetuleti). I am trying to match some
> user data. The are using ArcGis with the following definition:
>
>        PROJCS["HD_1972_Egyseges_Orszagos_Vetuleti",
>        GEOGCS["GCS_Hungarian_1972",
>        DATUM["D_Hungarian_1972",
>        SPHEROID["GRS_1967",6378160,298.247167427]],
>        PRIMEM["Greenwich",0],
>        UNIT["Degree",0.0174532925199432955]],
>        PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
>        PARAMETER["False_Easting",650000],
>        PARAMETER["False_Northing",200000],
>        PARAMETER["Scale_Factor",0.99993],
>        PARAMETER["Azimuth",90],
>        PARAMETER["Longitude_Of_Center",19.048571778],
>        PARAMETER["Latitude_Of_Center",47.14439372222],
>        UNIT["Meter",1]]
>
> I have set this up as:
>
>   from = pj_init_plus(
>                "+proj=latlong "
>                "+ellps=WGS84");
>
>   to = pj_init_plus(
>                "+proj=somerc "
>                "+lat_0=47.14439372222222 "
>                "+lon_0=19.04857177777778 "
>                "+k=0.99993 "
>                "+x_0=650000 "
>                "+y_0=200000 "
>                "+ellps=GRS67 "
>                "+units=m "
>                "+no_defs");
>
>     pj_transform(to, from, 1, 0, &le, &ln, &la);
>
> When I run some test data, I am always 80 some meters away. For example:
>
>   Ref:    Latitude:   47.7656110000 Longitude:   18.9532900000
>   Ref:    Easting:       642942.000 Northing:       269097.700
>   Calc:   Easting:       642857.473 Northing:       269067.872

By the way, a Hotine_Oblique_Mercator_Azimuth_Center is not the same as
somerc, the Swiss oblique Mercator projection. But, somerc is a better
approximation for EOV than a Hotine omerc.

Anyhow, with the datum transformation parameters from Mikael Rittri and a
coordinate frame rotation, I got with:
Ref:    Latitude:   47.7656110000 Longitude:   18.9532900000
a latitude, longitude of 47.765881, 18.954420
and with an accurate implementation of the EOV projection an (original axes)
x, y of:
269097.773, 642942.210
which is, let's say, spot on.
The approximation with somerc should give (its) x, y of:
642942.210, 269097.775


Oscar van Vlijmen


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

[Roger Oberholtzer]

On Tue, 2009-09-22 at 15:19 +0200, OvV_HN wrote:

> ----- Original Message -----
> From: "Roger Oberholtzer" <[hidden email]>
> To: <[hidden email]>
> Sent: Tuesday, September 22, 2009 11:02 AM
> Subject: [Proj] Problems with EOV for Hungary (1972 Egyseges
> OrszagosVetuleti)
>
>
> >I am trying to use proj (C API) to project WGS84 lat/longs in EOV for
> > Hungary (1972 Egyseges Orszagos Vetuleti). I am trying to match some
> > user data. The are using ArcGis with the following definition:
> >
> >        PROJCS["HD_1972_Egyseges_Orszagos_Vetuleti",
> >        GEOGCS["GCS_Hungarian_1972",
> >        DATUM["D_Hungarian_1972",
> >        SPHEROID["GRS_1967",6378160,298.247167427]],
> >        PRIMEM["Greenwich",0],
> >        UNIT["Degree",0.0174532925199432955]],
> >        PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
> >        PARAMETER["False_Easting",650000],
> >        PARAMETER["False_Northing",200000],
> >        PARAMETER["Scale_Factor",0.99993],
> >        PARAMETER["Azimuth",90],
> >        PARAMETER["Longitude_Of_Center",19.048571778],
> >        PARAMETER["Latitude_Of_Center",47.14439372222],
> >        UNIT["Meter",1]]
> >
> > I have set this up as:
> >
> >   from = pj_init_plus(
> >                "+proj=latlong "
> >                "+ellps=WGS84");
> >
> >   to = pj_init_plus(
> >                "+proj=somerc "
> >                "+lat_0=47.14439372222222 "
> >                "+lon_0=19.04857177777778 "
> >                "+k=0.99993 "
> >                "+x_0=650000 "
> >                "+y_0=200000 "
> >                "+ellps=GRS67 "
> >                "+units=m "
> >                "+no_defs");
> >
> >     pj_transform(to, from, 1, 0, &le, &ln, &la);
> >
> > When I run some test data, I am always 80 some meters away. For example:
> >
> >   Ref:    Latitude:   47.7656110000 Longitude:   18.9532900000
> >   Ref:    Easting:       642942.000 Northing:       269097.700
> >   Calc:   Easting:       642857.473 Northing:       269067.872
>
>
> By the way, a Hotine_Oblique_Mercator_Azimuth_Center is not the same as
> somerc, the Swiss oblique Mercator projection. But, somerc is a better
> approximation for EOV than a Hotine omerc.
>
> Anyhow, with the datum transformation parameters from Mikael Rittri and a
> coordinate frame rotation, I got with:
> Ref:    Latitude:   47.7656110000 Longitude:   18.9532900000
> a latitude, longitude of 47.765881, 18.954420
> and with an accurate implementation of the EOV projection an (original axes)
> x, y of:
> 269097.773, 642942.210
> which is, let's say, spot on.
> The approximation with somerc should give (its) x, y of:
> 642942.210, 269097.775

This looks promising. So, I have now tried:

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

   to = pj_init_plus(
       "+proj=somerc "
       "+lat_0=47.14439372222222 "
       "+lon_0=19.04857177777778 "
       "+k=0.99993 "
       "+x_0=650000 "
       "+y_0=200000 "
       "+ellps=GRS67 "
       "+units=m "
       "+towgs84=52.684,-71.194,-13.975,-0.312,-0.1063,-0.3729,1.0191 "
       "+no_defs");

     pj_transform(to, from, 1, 0, &le, &ln, &la);

And I see no difference at all. I use +towgs84 elsewhere, so I know it
works like this. But, as I said, I see no difference.

--
Roger Oberholtzer

OPQ Systems / Ramböll RST

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden

Office: Int [phone]
Mobile: Int [phone]

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

[Jean-Claude Repetto]

Roger Oberholtzer wrote :

>
>
> This looks promising. So, I have now tried:
>
>    from = pj_init_plus(
>        "+proj=latlong "
>        "+ellps=WGS84");
>
>    to = pj_init_plus(
>        "+proj=somerc "
>        "+lat_0=47.14439372222222 "
>        "+lon_0=19.04857177777778 "
>        "+k=0.99993 "
>        "+x_0=650000 "
>        "+y_0=200000 "
>        "+ellps=GRS67 "
>        "+units=m "
>        "+towgs84=52.684,-71.194,-13.975,-0.312,-0.1063,-0.3729,1.0191 "
>        "+no_defs");
>
>      pj_transform(to, from, 1, 0, &le, &ln, &la);
>
> And I see no difference at all. I use +towgs84 elsewhere, so I know it
> works like this. But, as I said, I see no difference.
>

You have forgotten to add the +towgs84 parameter to the first
projection. Alternatively, you can use "+datum=WGS84" instead of
"+ellps=WGS84 +towgs84=0,0,0,0,0,0,0".

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

[Roger Oberholtzer]

On Tue, 2009-09-22 at 17:25 +0200, Jean-Claude REPETTO wrote:

> You have forgotten to add the +towgs84 parameter to the first
> projection. Alternatively, you can use "+datum=WGS84" instead of
> "+ellps=WGS84 +towgs84=0,0,0,0,0,0,0".

This did indeed work. Thanks! I am now within 0.2 meters of the ArcGIS
calculations.

--
Roger Oberholtzer

OPQ Systems / Ramböll RST

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden

Office: Int [phone]
Mobile: Int [phone]

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