MITABCoordSys2SpatialRef not used adfDatumParm[0]-[6]

[warmerdam]

Dear Mr.. Warmerdam

I have used mitab_coordsys.cpp version 1.16 2000/10/16 21:44:50 and founded
that successfully tokenized and
stored in adfDatumParm[0]-[6] Bursa-Volf transformation parameters did not
used. If this not fixed today I can describe modifications I have made to
utilize Mif CoordSys parameters. I use maps on 1001 datum (Pulkovo 1942).

Best Regards
Sorry for my English.
Victor.

"Victor Osipkov" <vctos@email.com>

[warmerdam]

Frank,

Please note that at least this modification changes initial API description
{Before modification projdef->datum_params[6] was be in ?? units and after
modification PPM units supposed on input (according specification of
SetTOWGS84)} :

> I have inserted this conversion into pj_datum_set.c function: int
> pj_datum_set(paralist *pl, PJ *projdef) after operator
>         if( projdef->datum_params[3] != 0.0
>             || projdef->datum_params[4] != 0.0
>             || projdef->datum_params[5] != 0.0
>             || projdef->datum_params[6] != 0.0 ){
> as:
>         projdef->datum_params[6] = projdef->datum_params[6]/1000000.0+1;
> This is not the best place and you must desire where to place this code.

The problem is that this modification will really change PROJ API (datum
towgs84 definition procedure). I think one additional version uses angular
seconds as units of rotation parameters and ppm as scale parameter on input
of towgs required: there are many data sources and most of all uses angular
seconds for rotation parameters. Also as fool protection I recommend to
include into original function implementation error detection code for the
case of somebody will call this function using angular seconds instead of
radians (typical rotation in angular seconds is in 0.01-1 abs interval,
while in radians it will be something about 1e-8 - 1e-5). Idea of this
protection is: there many sources that not include units and not all users
understand that rotation is in angular seconds. What about me, I have lose 3
days because understand which units uses for rotation.

If you decide to change API I recommend also include another variant of 7
parameters Helmert transformation: Coordinate Frame Transformation (CFT).
CFT is a Bursa-Volf (BF) equivalent with backward rotation:
DxCFT = DxBF                    (datum_params[0])
DyCFT = DyBF                    (datum_params[1])
DzCFT = DzBF                    (datum_params[2])
RxCFT = -RxBF                   (datum_params[3])
RyCFT = -RyBF                   (datum_params[4])
RzCFT = -RzBF                   (datum_params[5])
MCFT = MBF                      (datum_params[6])

There are some information sources uses CFT instead BF as Helmert
transformation implementation and this simple function expand adaptability
of your advanced projects.

By the way: let me to thank you for sharing results of you job - it
dramatically decrease development time of tens geospatial projects!

Also I hope you will notify me if you found any errors and inaccuracy in my
changes.

Best Regards,
Victor Osipkov,
emailto:vctos@email.com
http://www.77.ru

[warmerdam]

1. In file mitab_coordsys.cpp
add lines at file begin:

        /*SEC_TO_RAD = Pi/180/3600*/
        #define SEC_TO_RAD 4.84813681109535993589914102357e-6

in OGRSpatialReference *MITABCoordSys2SpatialRef( const char * pszCoordSys )
after call SetGeogCS

        /* --------------------------------------------------------------------
*/
        /*      Set the GeogCS.
*/
        /* --------------------------------------------------------------------
*/
        poSR->SetGeogCS( pszGeogName, szDatumName, pszSpheroidName,
                     dfSemiMajor, dfInvFlattening,
                     pszPrimeM, dfPMLongToGreenwich,
                     SRS_UA_DEGREE,
                     atof(SRS_UA_DEGREE_CONV) );

insert call to SetTOWGS84

        poSR->SetTOWGS84( adfDatumParm[0], adfDatumParm[1], adfDatumParm[2],
                                                adfDatumParm[3]*SEC_TO_RAD,
                                                adfDatumParm[4]*SEC_TO_RAD,
                                                adfDatumParm[5]*SEC_TO_RAD,
                                                adfDatumParm[6] );

As I have understood, adfDatumParm values defined in MapInfoDatumInfo
asDatumInfoList[] is in angular seconds, while SetTOWGS84 require radians.
adfDatumParm[6] is steel in PPM while Bursa-Volf required scale factor M.
The conversion formula is

        M = 1+adfDatumParm[6]*1e-6

I have inserted this conversion into pj_datum_set.c function: int
pj_datum_set(paralist *pl, PJ *projdef) after operator
        if( projdef->datum_params[3] != 0.0
            || projdef->datum_params[4] != 0.0
            || projdef->datum_params[5] != 0.0
            || projdef->datum_params[6] != 0.0 ){
as:
        projdef->datum_params[6] = projdef->datum_params[6]/1000000.0+1;
This is not the best place and you must desire where to place this code.

2. Bursa-Volf calculation.
because I can not understand which units used in original implementation of
pj_geocentric_to_wgs84 I have replaced lines

            x_out = x[io] + defn->datum_params[0]
                + y[io] * defn->datum_params[5]
                + z[io] * defn->datum_params[4]
                + x[io] * defn->datum_params[6];

            y_out = y[io] + defn->datum_params[1]
                + x[io] * defn->datum_params[5]
                + z[io] * defn->datum_params[3]
                + y[io] * defn->datum_params[6];

            z_out = z[io] + defn->datum_params[2]
                + x[io] * defn->datum_params[4]
                + y[io] * defn->datum_params[3]
                + z[io] * defn->datum_params[6];

to lines:
            x_out = M_BF*(            x[io]  - Rz_BF*y[io] +
Ry_BF*z[io])+Dx_BF;
            y_out = M_BF*( Rz_BF*x[io] +            y[io]  -
Rx_BF*z[io])+Dy_BF;
            z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] +
z[io])+Dz_BF;

also I have added defines:
        #define Dx_BF defn->datum_params[0]
        #define Dy_BF defn->datum_params[1]
        #define Dz_BF defn->datum_params[2]
        #define Rx_BF defn->datum_params[3]
        #define Ry_BF defn->datum_params[4]
        #define Rz_BF defn->datum_params[5]
        #define M_BF  defn->datum_params[6]

We have:        Dx_BF, Dy_BF, Dz_BF in meters
                Rx_BF, Ry_BF, Ry_BF in radians
                M_BF as coefficient.

This is almost all.

For testing I have used the following code:
        OGRSpatialReference* srP = poLayer->GetSpatialRef();
        OGRSpatialReference wgsSR;
        OGRErr err = wgsSR.SetWellKnownGeogCS("WGS84");
        if (err != OGRERR_NONE){
                printf ("Unable to init WGS-84 %i\n", err);
        } else {
                transP = OGRCreateCoordinateTransformation( srP, &wgsSR);
        }

        if (transP){
                double x=5427843.85, y=7444438.88, z=0;
                printf("x=%.1f, y=%.1f, z=%.1f -> ", x, y, z);
                transP->Transform(1, &x, &y, &z);
                printf("x=%.6f, y=%.6f, z=%.6f\n", x, y, z);
        }


I nave used mif file with the following header:

Version 300
Charset "WindowsCyrillic"
Delimiter ","
CoordSys Earth Projection 8, 1001, "m", 56.03333, 67.05, 1, 5400000, 7429000
Bounds (-2849281.53901, -10013196.8816) (13649281.539, 9991078.11391)
Columns 1
  autocad_elevation Decimal(9, 7)
Data

MapInfo convert (5427843.85, 7444438.88, 0) as WGS84:(56.67467086,
67.18761636, ??)
test convert (5427843.85, 7444438.88, 0) as WGS84:(56.674488,
67.187530, -4.1)
The position calculation difference is 12m.
I suppose this error produced because conversion parameters is not the same
as used in MapInfo and we have converted points near Polar circle. What are
you think about it?

Best Regards,
Victor

[warmerdam]

Essentially all changes applied, though I have decided that the values
passed to SetTOWGS84 should have angles in arc seconds, and scaling
values in ppm.  PROJ.4, MITAB and OGRSpatialReference updated accordingly.

[warmerdam]

Victor,

I have finally returned to the TOWGS84 issue, and I think I will need
to have the rotational parameters in arc seconds and the scaling factor
in "parts per million" to adhere to the OGC specifications.  This would
correspond to the EPSG transformation type called "Position Vector 7-param
transformation" (COORD_TRF_METHOD_CODE=9606).

I will add some of your suggested argument checking.

Folks with information in Coordinate Frame Rotation will be expected to
reverse the signs on the rotation parameters before using with
OGRSpatialReference and PROJ.4.

I really appreciate your comprehensive explanation of the problems and
issues.

Best regards,