incorrect data for EPSG:3120 in FWTools2.4.2 and 2.4.6

[dkrysian]

Hello, I'm sorry for my english.... Problem occurs while converting data between EPSG:3120 and EPSG:4326 - output geometry is moved ~100m from correct position. I must use correct datafiles from FWTools2.2.8.

And another problem: when converting data between EPSG:3120 and EPSG:2178 ogr2ogr always return incorrect position, but when converting EPSG:3120 --> EPSG:4326 and then EPSG:4326 --> EPSG:2178 position is correct! anybody knows why?

[Even Rouault]

I CC Frank as being the expert for how the .CSV files are generated from the EPSG database and Andrey for his particular interest in Pulkovo ;-)

My findings are the following one :

EPSG:3120 is PROJCSPulkovo 1942(58) / Poland zone I which relies on GCSPulkovo 1942(58) of code 4179

data/gcs.csv in GDAL 1.7.0 contains : 4179,"Pulkovo 1942(58)",6179,"Pulkovo 1942(58)",6179,9122,7024,8901,1,0,6422

whereas the same file in FWTools 2.2.8 contains : 4179,"Pulkovo 1942(58)",6179,"Pulkovo 1942/58",6179,9122,7024,8901,1,0,9606,33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84

So there is no more the datum shift parameters in newer GDAL data files, which explains the offset of the results. If this is a feature or a bug, I can't say. Generally this means that there was some ambiguity in the original EPSG database concerning the datum shift parameters, and thus they have been dropped when converting to GDAL .CSV files.

Maybe this could be dealt as an "override", as done for #3176, but this is far beyond my level of expertise to say if it is appropriate or not.

In FWTools 2.2.8, testepsg EPSG:3120 expands to the following PROJ.4 string : "+proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 +units=m +no_defs". So if you use that string for specifying the SRS, newer GDAL should give identical results to FWTools 2.2.8, like :

ogr2ogr -s_srs "+proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 +units=m +no_defs" -t_srs EPSG:4326 out_4326.shp in_3120.shp

[warmerdam]

Interesting - my earlier answer on the (fwtools?) mailing list may have been misleading.

This issue - the transform disappearing in a newer version - can happen if EPSG goes from publishing one transformation to publishing multiple transformations. When this happens, I cannot choose which is best so instead I choose none of them.

I will revisit the generation techniques at some point in the future.

[dron]

Replying to warmerdam:

This issue - the transform disappearing in a newer version - can happen if EPSG goes from publishing one transformation to publishing multiple transformations. When this happens, I cannot choose which is best so instead I choose none of them.

I'm afraid that we can't make such a decision without the end user. Only user knows what is the best transformation for his data. Or she doesn't know it and will try all the options available. I think there should be a way to choose among multiple transformations (just like ESRI software does).

[dkrysian]

Thanks for reply, this problem affects to 5 Poland zones (3120,2172,2173,2174,2175) on Krassowski ellipsoid - for this zones transformation data in FWTools2.2.8 toWgs84 was determined and was correct. I think, every EPSG ID should have determined one transformation to WGS84.

[warmerdam]

Andrey,

Can you clarify how ESRI software does this? I don't believe I have ever attempted to do a coordinate transformation with ESRI software. Does it pop up a dialog with all the known datum shift methods for a given datum?

I'm not sure how we would accomplish this in the context of GDAL. We could attempt to capture all the transformations from EPSG, mark one as the default, and offer a mechanism to list all the shifts available with a given datum and to alter the default though this would be fairly involved.

[dron]

Replying to warmerdam:

Can you clarify how ESRI software does this? I don't believe I have ever attempted to do a coordinate transformation with ESRI software. Does it pop up a dialog with all the known datum shift methods for a given datum?

Exactly. When you load dataset with the datum different from the one specified for your workspace ArcGIS shows a dialog with a list of available transformations including an option to omit the transformation. Though I can't call their implementation excellent because it lacks one important feature: description of the particular transformation. To pick the right one you should read the coefficients for each one and decide which is the best for you using some external knowledge. If we will do that in GDAL/OGR we should keep the descriptions from the EPSG database in our tables and make them available for user to read when choosing the transformation. In majority of cases there will be enough information to do the right decision.

I'm not sure how we would accomplish this in the context of GDAL. We could attempt to capture all the transformations from EPSG, mark one as the default, and offer a mechanism to list all the shifts available with a given datum and to alter the default though this would be fairly involved.

Agreed I am not sure how to implement it in the best way. Maybe OGRCreateCoordinateTransformation() should fail if there are multiple datum shifts available for supplied datums? It can take the third parameter with the datum shift of a choice. Also the datum shift can be set by a specific method of OGRCoordinateTransformation. Whatever the case this should be implemented in some way in PROJ.4 as a bottom line.

[warmerdam]

As part of the "big datum shift upgrade" in libgeotiff, we can now record all the datum shifts available fora GCS in the datum_shift.csv file. The one marked as preferred is used in the gcs.csv file. Currently the datum_shift_pref.csv file is used to set operation 1267 as the preferred solution for Pulkovo 1942 with these parameters:

        TOWGS84[33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84],

Does this seem reasonable?

While the datum shift options are available in datum_shift.csv there is currently no user interface options available to show or select one.

[dron]

Replying to warmerdam:

As part of the "big datum shift upgrade" in libgeotiff, we can now record all the datum shifts available fora GCS in the datum_shift.csv file. The one marked as preferred is used in the gcs.csv file. Currently the datum_shift_pref.csv file is used to set operation 1267 as the preferred solution for Pulkovo 1942 with these parameters: TOWGS84[33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84]

Does this seem reasonable?

Yes, it does. 1267 is fine. The only note from my side that it's parameters are

TOWGS84[23.92,-141.27,-80.9,-0,0.35,0.82,-0.12]

and not the ones you've listed above.

Thank you!

[warmerdam]

Sorry, I got mixed up between Pulkovo_1942_58 and Pulkovo_1942. The parameters I quoted are for GCS EPSG:4179 - Pulkovo_1942_58 which is used by EPSG:3357 which this ticket is about. The value you quote is selected for Pulkovo_1942 (EPSG:4284).

I think we are good ... let me know if I'm off base.

WKT[EPSG:4179] =
GEOGCS["Pulkovo 1942(58)",
    DATUM["Pulkovo_1942_58",
        SPHEROID["Krassowsky 1940",6378245,298.3,
            AUTHORITY["EPSG","7024"]],
        TOWGS84[33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84],
        AUTHORITY["EPSG","6179"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4179"]]

WKT[EPSG:4284] =
GEOGCS["Pulkovo 1942",
    DATUM["Pulkovo_1942",
        SPHEROID["Krassowsky 1940",6378245,298.3,
            AUTHORITY["EPSG","7024"]],
        TOWGS84[23.92,-141.27,-80.9,-0,0.35,0.82,-0.12],
        AUTHORITY["EPSG","6284"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4284"]]