[Qgis-user] Re: Annoying CRS-Problem with EPSG 31468 / 2167: Conclusions?

Sun Apr 1 23:27:58 PDT 2012

Am 30.03.2012 19:44, schrieb Etienne Tourigny:
> On Fri, Mar 30, 2012 at 6:23 AM,<bernhard.stroebl at jena.de>  wrote:
>> Etienne,
>>
>> I see this stuff from a user's point of view. I am currently not working
>> with different projections but my projects are EPSG:31468 and my data, too.
>> So this is my perspective. I am neither a geodesist nor did I dig into how
>> different software handles srs or prj files. All that I know is that QGIS
>> (1.7.4) imports EPSG:31468 shape files as EPSG:2167 because it does neither
>> evaluate the DATUM["D_Deutsches_Hauptdreiecksnetz" nor the
>> SPHEROID["Bessel_1841" infos given.
>> Trunk creates a user-defined SRS and does not apply EPSG:31468 either.
>> I would prefer QGIS to recognize the shape file as being EPSG:31468.
>
> Strictly speaking it is not exactly EPSG:31468 (that would require the
> AUTHORITY node, which .prj files do not have), but the PROJ.4 string
> should match that of EPSG:31468, at least.

yup

> That is the problem with .prj files, it`s not easy to relate them to EPSG.

then I would like to come back to my former suggestion: analyze the prj 
file and find any matching EPSG-code
one EPSG -> OK but tell user
several -> let user choose (and probably remember choice)
none -> ask user
>
>>
>>
>> Am 30.03.2012 03:23, schrieb Etienne Tourigny:
>> [snip]
>>
>>
>>>>
>>>>>
>>>>> You can also play with the command line and see what comes up.
>>>>>
>>>>> GDAL_FIX_ESRI_WKT=GEOGCS gdalsrsinfo file.prj
>>>>
>>>>
>>>>
>>>> PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
>>>> +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m
>>>> +no_defs '
>>>>
>>>>> GDAL_FIX_ESRI_WKT=TOWGS84 gdalsrsinfo file.prj
>>>>
>>>>
>>>>
>>>> PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
>>>> +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m
>>>> +no_defs '
>>>>
>>>>> GDAL_FIX_ESRI_WKT=NO gdalsrsinfo file.prj
>>>>
>>>>
>>>>
>>>> PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
>>>> +datum=potsdam +units=m +no_defs '
>>>>
>>>> This would be my preferred parameters as datum=potsdam can be different
>>>> towgs84 parameters in different locations
>>>
>>>
>>> How would this be decided? By the user or some magic way?
>>
>>
>> The user can change it with the plugin "Coordinate Systems Updater". You
>> see: if the datum gets translated into a towgs84 parameter and the
>> projection is matched based on this then if a different user-defined towgs84
>> paramter exists the right projection is not found.
>> according to http://mapref.org/GeodeticReferenceSystemsDE.html#Zweig733 the
>> towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7
>> is applicable for the whole country of Germany
>> but local towgs84 parameters may differ
>>
>
> Then it makes sense to have those +towgs84 parameters as the default,
> if you leave them out then you get an incorrect datum shift.

agreed
>
>>
>>>
>>> This sort of contradicts your statement that it should be recognized
>>> as EPSG:31468 because that code is NOT potsdam datum (in GDAL/OGR
>>> anyway)
>>
>>
>> AFAIK the DHDN is the realisation of the Potsdam Datum see e.g. here:
>> http://georepository.com/datum_6314/Deutsches-Hauptdreiecksnetz.html
>> so in the prj file it is called DATUM["Deutsches_Hauptdreiecksnetz" whereas
>> Proj4 says +datum=potsdam ; both mean the same
>
> OK
>
> Using proj-4.7.1 and gdal-1.9 the +datum=potsdam gives different
> towgs84 parameters than EPSG:4314. This has been fixed in proj-4.8.
> Could this contribute to the problem?

I would not think so, as it is just another towgs84 parameter set which 
gives a different position when making a 7-parameter transformation
>
> $ gdalsrsinfo '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
> +datum=potsdam +units=m +no_defs '
>
> PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
> +ellps=bessel +towgs84=606,23,413,0,0,0,0 +units=m +no_defs '
>
> $ gdalsrsinfo EPSG:4314
>
> PROJ.4 : '+proj=longlat +ellps=bessel
> +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +no_defs '
>
>>
>>>
>>>>
>>>> However the problem persists that QGIS does not recognize the projection
>>>> as
>>>> EPSG:31468 no matter if I set it to GEOGCS, TOWGS84 or NO!
>>>> Neiter does it detect the correct projction from the prj file QGIS itself
>>>> created.
>>>
>>>
>>> The GEOGCS option to GDAL_FIX_ESRI_WKT fixes the GEOGCS authority
>>> node, but not the root node which is a PROJCS (see example below)!
>>> Trying to match all possible projections can be more challenging than
>>> matching the possible geogcs definitions. It is possible, but I didn't
>>> get there yet.
>>>
>>> Just to make sure I understand - is the expected behaviour in QGis
>>> (from you POV) to see the file as EPSG:31468 with the TOWGS84
>>> parameters seen above? Or the more generic 'potsdam' datum?
>>
>>
>> Maybe it would be better to see the towgs84 parameters (as defined by QGIS
>> default or the user via the plugin) but
>> 1) it must be imported as EPSG:31468 (and not as user-defined)
>> 2) if saving to shapefile only the DATUM["Deutsches_Hauptdreiecksnetz" must
>> be written (as QGIS does) and IMHO not the
>> TOWGS84[598.1,73.7,418.2,0.202,0.045,-2.455,6.7]
>
> I agree. What is written to .prj file by QGis when you choose EPSG:31468 ?

1.7.4: 
PROJCS["DHDN_3_degree_Gauss_Kruger_zone_4",GEOGCS["GCS_DHDN",DATUM["D_Deutsches_Hauptdreiecksnetz",
SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",12],
PARAMETER["scale_factor",1],PARAMETER["false_easting",4500000],PARAMETER["false_northing",0],
UNIT["Meter",1]]

which gets imported as EPSG:2167 when deleting the qpj file.

>
>>
>>
>>>
>>> ===>
>>> Perhaps QGis could be modified to scan existing proj4 definitions for
>>> a match with the layer's proj4 string as Berhard suggests.
>>> That should probably work in this case (because all that is missing is
>>> the authority nodes - which is probably the comparison inside Qgis) -
>>> am I right Jurgen?
>>>
>>>
>>> Example of working GEOGCS detection:
>>>
>>> $  GDAL_FIX_ESRI_WKT=GEOGCS gdalsrsinfo  tmp1.prj
>>>
>>> PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
>>> +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7
>>> +units=m +no_defs '
>>>
>>> OGC WKT :
>>> PROJCS["Transverse_Mercator",
>>>      GEOGCS["DHDN",
>>>          DATUM["Deutsches_Hauptdreiecksnetz",
>>>              SPHEROID["Bessel 1841",6377397.155,299.1528128,
>>>                  AUTHORITY["EPSG","7004"]],
>>>              TOWGS84[598.1,73.7,418.2,0.202,0.045,-2.455,6.7],
>>>              AUTHORITY["EPSG","6314"]],
>>>          PRIMEM["Greenwich",0,
>>>              AUTHORITY["EPSG","8901"]],
>>>          UNIT["degree",0.0174532925199433,
>>>              AUTHORITY["EPSG","9122"]],
>>>          AUTHORITY["EPSG","4314"]],
>>>      PROJECTION["Transverse_Mercator"],
>>>      PARAMETER["latitude_of_origin",0],
>>>      PARAMETER["central_meridian",12],
>>>      PARAMETER["scale_factor",1],
>>>      PARAMETER["false_easting",4500000],
>>>      PARAMETER["false_northing",0],
>>>      UNIT["Meter",1]]
>>
>>
>> this is intersting as the definition has a datum and a towgs84 parameter
>> AFAIK the towgs84 values are translating the datum for the use in a
>> 7-parameter datum transformation. This becomes relevant when my project has
>> a srs using a WGS84 ellipsoid and I need otf transformation of data defined
>> in EPSG:31468
>> I doubt it is useful to store a DATUM _and_ a TOWGS84 node, only the datum
>> is part of the srs and thus part of the data whereas the question how to
>> transform to wgs84 is more the responsibility of the software that performs
>> the transformation.
>
> in GDAL/OGR if you do not have TOWGS84 then the projection is
> effectively treated as TOWGS84=0,0,0,0,0,0,0 so you do need a TOWGS84
> node.  In proj.4 that is a different story (I think).
> So internally, you do need to have TOWGS84 node in addition to DATUM.

OK

> See the CT-1.0 specs.
>
>
>>
>>
>>
>>>
>>> and compare with EPSG:31468, the GEOGCS node has correct authority (4314)
>>>
>>> $ gdalsrsinfo EPSG:31468
>>>
>>> PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
>>> +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7
>>> +units=m +no_defs '
>>>
>>> OGC WKT :
>>> PROJCS["DHDN / 3-degree Gauss-Kruger zone 4",
>>>      GEOGCS["DHDN",
>>>          DATUM["Deutsches_Hauptdreiecksnetz",
>>>              SPHEROID["Bessel 1841",6377397.155,299.1528128,
>>>                  AUTHORITY["EPSG","7004"]],
>>>              TOWGS84[598.1,73.7,418.2,0.202,0.045,-2.455,6.7],
>>>              AUTHORITY["EPSG","6314"]],
>>>          PRIMEM["Greenwich",0,
>>>              AUTHORITY["EPSG","8901"]],
>>>          UNIT["degree",0.0174532925199433,
>>>              AUTHORITY["EPSG","9122"]],
>>>          AUTHORITY["EPSG","4314"]],
>>>      PROJECTION["Transverse_Mercator"],
>>>      PARAMETER["latitude_of_origin",0],
>>>      PARAMETER["central_meridian",12],
>>>      PARAMETER["scale_factor",1],
>>>      PARAMETER["false_easting",4500000],
>>>      PARAMETER["false_northing",0],
>>>      UNIT["metre",1,
>>>          AUTHORITY["EPSG","9001"]],
>>>      AXIS["X",NORTH],
>>>      AXIS["Y",EAST],
>>>      AUTHORITY["EPSG","31468"]]
>>>
>>> Etienne
>>>
>>
>> [more snip]
>>
>>
>> How does QGIS work? What happens when it is confronted with a prj file? How
>> does the matching work if no authority is given for the PROJCS?
>
> I have no idea here, hoping some of the qgis devs would comment.
> Perhaps sending a mail to the qgis-dev list is more proper.
>
>> With EPSG:31468 we have the problem that the srs string uses a different
>> name for the datum than Proj4 and that this datum can obtain different
>> towgs84 paramters for different parts of the world.
>
> I might add this problem is not unique for your case.  In fact many
> datums have different datum shifts depending on location. GDAL/OGR
> uses the "preferred" datum shift values.

which is fine

> I think that this is a different issue, and is up to the user to
> select which datum is needed on a per-case basis (or change the
> default values as you suggested).

agreed, I was just trying to sum up open questions in relation to QGIS

>
> What needs to be resolved is the .prj<->WKT/PROJ.4 relations and
> correct preffered datum and EPSG code (or Qgis representation of
> them).

agreed!

>
>>
>> Bernhard
>>
>>
>>
>> ________ Information from NOD32 ________
>> This message was checked by NOD32 Antivirus System for Linux Mail Server.
>> http://www.nod32.com

________ Information from NOD32 ________
This message was checked by NOD32 Antivirus System for Linux Mail Server.
http://www.nod32.com