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

Etienne Tourigny etourigny.dev at gmail.com
Thu Mar 29 18:23:05 PDT 2012 

On Mon, Mar 26, 2012 at 3:42 AM,  <bernhard.stroebl at jena.de> wrote:
> Am 23.03.2012 13:32, schrieb Etienne Tourigny:
>
>> On Fri, Mar 23, 2012 at 5:26 AM,<bernhard.stroebl at jena.de>  wrote:
>>>
>>> I just tested with current trunk (compiled on OpenSUSE 64bit):
>>> When I load
>>>
>>> PROJCS["Transverse_Mercator",GEOGCS["GCS_bessel",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 is EPSG:31468 (proj4 params in QGIS: +proj=tmerc +lat_0=0 +lon_0=12
>>> +k=1 +x_0=4500000 +y_0=0 +ellps=bessel
>>> +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs)
>>> QGIS creates a user-defined srs:
>>> +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
>>>
>>> although there is a difference in towgs84 params the layer matches a
>>> corresponding PostGIS layer defined with EPSG:31468 (I am aware there
>>> exist
>>> different towgs84-parameter sets for this datum)
>>>
>>> Same happens when saving as shape file from QGIS deleting the qpj file
>>> and
>>> then reloading into QGIS, so QGIS is not able to correctly interpret its
>>> own
>>> prj file!
>>>
>>> The only way I found to get the shape file loaded as EPSG:31468 is
>>> appending>>AUTHORITY["EPSG","31468"]<<  to the prj file
>>>
>>> PostGIS' spatial_ref_sys table contains the prj string (field srtext) and
>>> the proj4 params (field proj4text) for the projection
>>> srtext: PROJCS["DHDN / 3-degree Gauss-Kruger zone
>>> 4",GEOGCS["DHDN",DATUM["Deutsches_Hauptdreiecksnetz",SPHEROID["Bessel
>>>
>>> 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6314"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4314"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",12],PARAMETER["scale_factor",1],PARAMETER["false_easting",4500000],PARAMETER["false_northing",0],AUTHORITY["EPSG","31468"],AXIS["Y",EAST],AXIS["X",NORTH]]
>>> proj4text: +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
>>> +ellps=bessel +datum=potsdam +units=m +no_defs
>>>
>>> remarks:
>>> 1) srtext contains>>AUTHORITY["EPSG","31468"]<<
>>> 2) proj4 says>>+datum=potsdam<<  whereas srtext says
>>>>>
>>>>> DATUM["Deutsches_Hauptdreiecksnetz"<<
>>>
>>>
>>> What could help?
>>> Somehow QGIS needs to use _all_ information given in the prj file and
>>> find a
>>> matching srs(proj4 parameters) for that. If there are several matching
>>> the
>>> user should choose which to apply, if there is none matching create a
>>> user-defined srs with the parameters given but tell the user. Maybe the
>>> datum parameter should be added to the proj4 definition of the
>>> projection. I
>>> am aware that a +datum parameter is only a short for a +towgs84 parameter
>>> but the datum is the same for the projection whereas towgs84 may differ
>>> locally.
>>>
>>> Maybe in at most 5 years' time we won't need GAUSS-KRUEGER any more as we
>>> are supposed to be switched to ETRS89 by then ;-)
>>>
>>> Bernhard
>>>
>>> Am 22.03.2012 19:23, schrieb Etienne Tourigny:
>>>
>>>> The following only applies to builds using master (and upcoming 1.8)
>>>> and gdal-1.9
>>>>
>>>> Jeff pushed in master the fix which add missing TOWGS84 parameters (by
>>>> using GDAL_FIX_ESRI_WKT=TOWGS84).  It is possible that using
>>>> GDAL_FIX_ESRI_WKT=GEOGCS would resolve this issue.  Perhaps we couls
>>>> add an option for this, and set the default to GEOGCS?
>>>
>>>
>>>
>>> how can I use this?
>>
>>
>> Berhard, this is already in Master and 1.8.  Can you try setting
>> GDAL_FIX_ESRI_WKT=GEOGCS in qgsogrprovider.cpp?
>
>
> Etienne,
>
> thank you for your instructions
> I replaced WGS84 in qgsogrprovider.cpp with GEOGCS (QGIS master)
> result stays the same: User-defined srs with towgs84 params is applied
>
>
>>
>> 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?

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)

>
> 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?

===>
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]]

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

>
> Bernhard
>
>
>>
>> although in this case I am afraid that gdal/ogr sees it as
>> Deutsches_Hauptdreiecksnetz
>>
>> all the best,
>> Etienne
>
>
>
>
>
>>>
>>>>
>>>> see bug http://hub.qgis.org/issues/5142
>>>>
>>>> Is there a proper fix for this in 1.7.4?
>>>>
>>>> Etienne
>>>>
>>>>
>>>> On Thu, Mar 22, 2012 at 12:35 PM,<bernhard.stroebl at jena.de>    wrote:
>>>>>
>>>>>
>>>>>
>>>>> Am 22.03.2012 16:05, schrieb Bernd Vogelgesang:
>>>>>
>>>>>> So, the problem was now described thoroughly, but what are the
>>>>>> conclusions now?
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> There is already a ticket open http://hub.qgis.org/issues/4977
>>>>> although I have no idea if the issue has been solved or not
>>>>>
>>>>>>
>>>>>> I'm going to spread QGIS in my region in the next months, and it's
>>>>>> quite
>>>>>> hard to explain to the people, who to 99% work with EPSG 31468-data,
>>>>>> that they can't ever rely on any data they produce and that they have
>>>>>> to
>>>>>> double-check each layer after creation for the right CRS, even if they
>>>>>> set it correctly in each setting available?!?
>>>>>>
>>>>>> I'm not too familiar with all those packages involved handling the
>>>>>> CRSs,
>>>>>> so can someone please point me the way who to address with this
>>>>>> problem?
>>>>>>
>>>>>> For novice QGIS-users, this "bug" is a real exclusion criterion for
>>>>>> QGIS
>>>>>> for they will only create crap with it!
>>>>>>
>>>>>> Again my findings:
>>>>>> Options settings: new layers and projects are created with EPSG 31468,
>>>>>> promting for CRS ist set.
>>>>>>
>>>>>> - EPSG 31468-Layers with a prj-file from ArcGIS get loaded as 2167
>>>>>> without promting for the CRS
>>>>>> - EPSG 31468-Layers loaded into a EPSG 31468-Region in GRASS are
>>>>>> handed
>>>>>> back to QGIS as ... 2167
>>>>>> - New Layers created in some plugins (like some ftools-plugins), added
>>>>>> to the TOC, are in ... 2167 .. or even 3397
>>>>>>
>>>>>>
>>>>>> What can i do? (And i really wonder what all those other german users
>>>>>> do
>>>>>> and why the don't complain...)
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> For me the base of the problem seems to be that the prj file is not
>>>>> matched
>>>>> to the right EPSG srs. As soon as you tell the layer it is EPSG:31468
>>>>> it
>>>>> is
>>>>> positioned correctly. So the problem arises only with shapefiles; at my
>>>>> work
>>>>> we use PostGIS for most of our data (that is of course no help for you
>>>>> but
>>>>> explains why I am not "complaining")
>>>>>
>>>>> Bernhard
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> Thanx for advice
>>>>>> Bernd
>>>>>>
>>>>>>
>>>>>> (QGIS 1.7.4, osgeo4w advanced install, Win7 64bit)
>>>>>>
>>>>>>
>>>>>>
>>>
>>>
>>> ________ Information from NOD32 ________
>>> This message was checked by NOD32 Antivirus System for Linux Mail Server.
>>> http://www.nod32.com
>
>
> --
> Bernhard Ströbl
> Anwendungsbetreuer GIS
>
> Kommunale Immobilien Jena
> Am Anger 26
> 07743 Jena
>
> Tel.: 03641 49- 5190
> E-Mail: bernhard.stroebl at jena.de
> Internet: www.kij.de
>
>
>
> Kommunale Immobilien Jena
> Eigenbetrieb der Stadt Jena
> Werkleiter: Thomas Dirkes
>
>
>
> ________ Information from NOD32 ________
> This message was checked by NOD32 Antivirus System for Linux Mail Server.
> http://www.nod32.com

More information about the Qgis-user mailing list