Opened 23 years ago
Last modified 23 years ago
#41 closed defect (fixed)
MITABCoordSys2SpatialRef not used adfDatumParm[0]-[6]
Reported by: | warmerdam | Owned by: | warmerdam |
---|---|---|---|
Priority: | high | Milestone: | |
Component: | default | Version: | unspecified |
Severity: | normal | Keywords: | |
Cc: |
Description
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>
Change History (4)
comment:2 by , 23 years ago
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
comment:3 by , 23 years ago
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.
comment:4 by , 23 years ago
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,
Note:
See TracTickets
for help on using tickets.