[GRASS5] 5.7: Datum transformation - DHDN?

GRASS GIS GRASS GIS developer
1

Hello Paul
(cc grass5)

I'm just having import fun with projections and would like
to hear your opinion:

Trying to import the German FRIDA free vector data extended by a .prj
file which looks like this:

cat frida-1.0-shp-joined/proj.prj
PROJCS["Transverse Mercator",GEOGCS["bessel",DATUM["D_Deutsche_Hauptdreiecksnetz",SPHEROID["bessel",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",1],PARAMETER["false_easting",3500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

I face the following problem:

v.in.ogr dsn=frida-1.0-shp-joined/strassen-joined.shp output=strassen location=frida
A datum name potsdam (Deutsches_Hauptdreiecksnetz) was specified without transformation parameters.

Now select Datum Transformation Parameters
Enter 'list' to see the list of available Parameter sets
Enter the corresponding number, or <RETURN> to cancel request

Checking the name, reports:
grep -i Hauptdreiec ~/grass57/lib/gis/*
~/grass57/lib/gis/datum.table:potsdam "Deutsches_Hauptdreiecksnetz" bessel dx=606.0 dy=23.0 dz=413.0

Since ESRI uses the (gramatically wrong) 'D_Deutsche_Hauptdreiecksnetz', GRASS fails to
find the values. So far no surprise.

To overcome the problem I would like to add 'D_Deutsche_Hauptdreiecksnetz'
in lib/gis/datum.table. But I assume that I should add the 7 parms version
which I could extract from the EPSG SQL database
(should be identical to 'DHDN' = Deutsches HauptDreiecksNetz):

I created a VIEW 'datum_parameters_tidied' like this:

SELECT o.coord_op_name, m.coord_op_method_name, p.parameter_name, pu.sort_order, pv.parameter_value, u.unit_of_meas_name, a.area_name FROM ((((((epsg_coordoperation o JOIN epsg_coordoperationmethod m ON ((m.coord_op_method_code = o.coord_op_method_code))) JOIN epsg_coordoperationparamvalue pv ON (((pv.coord_op_method_code = o.coord_op_method_code) AND (pv.coord_op_code = o.coord_op_code)))) JOIN epsg_coordoperationparamusage pu ON (((pu.coord_op_method_code = o.coord_op_method_code) AND (pu.parameter_code = pv.parameter_code)))) JOIN epsg_coordoperationparam p ON ((p.parameter_code = pv.parameter_code))) JOIN epsg_unitofmeasure u ON ((u.uom_code = pv.uom_code))) JOIN epsg_area a ON ((a.area_code = o.area_of_use_code))) WHERE (((o.coord_op_type = 'transformation'::character varying) AND (o.coord_op_name ~~ '%WGS 84%'::text)) AND (o.deprecated = 0)) ORDER BY o.coord_op_name;

which then reports:
select * from datum_parameters_tidied where coord_op_name ~ 'DHDN';

[ see attached file dhdn_datum_parameters.csv ]

Question: which one to pick? I assume the
'Position Vector 7-param. transformation'...

Recommendations are welcome.

Thanks
Markus

(attachments)

dhdn_datum_parameters.csv (1.9 KB)

2

Hello Markus

On Wed, 6 Oct 2004, Markus Neteler wrote:

Hello Paul
(cc grass5)

I'm just having import fun with projections and would like
to hear your opinion:

Trying to import the German FRIDA free vector data extended by a .prj
file which looks like this:

cat frida-1.0-shp-joined/proj.prj
PROJCS["Transverse Mercator",GEOGCS["bessel",DATUM["D_Deutsche_Hauptdreiecksnetz",SPHEROID["bessel",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",1],PARAMETER["false_easting",3500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

I face the following problem:

v.in.ogr dsn=frida-1.0-shp-joined/strassen-joined.shp output=strassen location=frida
A datum name potsdam (Deutsches_Hauptdreiecksnetz) was specified without transformation parameters.

Now select Datum Transformation Parameters
Enter 'list' to see the list of available Parameter sets
Enter the corresponding number, or <RETURN> to cancel request

^^^^^^^^^

If you type list here it should give you a choice of the GRASS parameters for the DHDN/Potsdam datum, I think.

Checking the name, reports:
grep -i Hauptdreiec ~/grass57/lib/gis/*
~/grass57/lib/gis/datum.table:potsdam "Deutsches_Hauptdreiecksnetz" bessel dx=606.0 dy=23.0 dz=413.0

Since ESRI uses the (gramatically wrong) 'D_Deutsche_Hauptdreiecksnetz', GRASS fails to
find the values. So far no surprise.

But note that the warning message from GRASS used the correct format. I believe this is due to some good work by Frank in the OSRMorphFromESRI() from function which is called by GRASS to fix things up before importing the projection information, so I don't think there is a problem here really.

To overcome the problem I would like to add 'D_Deutsche_Hauptdreiecksnetz'
in lib/gis/datum.table.

The 'correct' way to do this (although it's really a bit of a hack) is to add two lines to the static const char *papszDatumEquiv in lib/proj/convert.c and recompile. The first one should be the 'wrong' value that has occured in your data, and the second one the value in the
GRASS datum.table. I did this for the New Zealand 2000 datum that Hamish was having some trouble with. Probably the even more correct way to fix it is to ask Frank to add some code to OGR to fix it.

But as I said I don't think it's necessary here---please check.

Paul

3

Hello Paul,

On Wed, Oct 06, 2004 at 03:18:37PM +0100, Paul Kelly wrote:

Hello Markus

On Wed, 6 Oct 2004, Markus Neteler wrote:

>Hello Paul
>(cc grass5)
>
>I'm just having import fun with projections and would like
>to hear your opinion:
>
>Trying to import the German FRIDA free vector data extended by a .prj
>file which looks like this:
>
>cat frida-1.0-shp-joined/proj.prj
>PROJCS["Transverse
>Mercator",GEOGCS["bessel",DATUM["D_Deutsche_Hauptdreiecksnetz",SPHEROID["bessel",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",1],PARAMETER["false_easting",3500000],PARAMETER["false_northing",0],UNIT["Meter",1]]
>
>I face the following problem:
>
>v.in.ogr dsn=frida-1.0-shp-joined/strassen-joined.shp output=strassen
>location=frida
>A datum name potsdam (Deutsches_Hauptdreiecksnetz) was specified without
>transformation parameters.
>
>Now select Datum Transformation Parameters
>Enter 'list' to see the list of available Parameter sets
>Enter the corresponding number, or <RETURN> to cancel request
>>
^^^^^^^^^

If you type list here it should give you a choice of the GRASS parameters
for the DHDN/Potsdam datum, I think.

Yes, right:

list

Number Details
---
1 Used in Germany (Sachsen)
        (PROJ.4 Params towgs84=612.4,77.0,440.2,-0.054,0.057,-2.797,2.55)
        Accuracy <1m
---
2 Used in Germany (Thüringen)
        (PROJ.4 Params towgs84=599.4,72.4,419.2,-0.063,-0.022,-2.723,6.46)
        Accuracy <1m
---
3 Used in Germany (West - North - 52d20'N to 55d00'N)
        (PROJ.4 Params towgs84=590.5,69.5,411.6,-0.796,-0.052,-3.601,8.30)
        Accuracy <1m
---
4 Used in Germany (West - Middle - 50d20'N to 52d20'N)
        (PROJ.4 Params towgs84=584.8,67.0,400.3,0.105,0.013,-2.378,10.29)
        Accuracy <1m
---
5 Used in Germany (West - South - 47d00N to 50d20'N)
        (PROJ.4 Params towgs84=597.1,71.4,412.1,0.894,0.068,-1.563,7.58)
        Accuracy <1m
---
6 Used in Germany (Whole Country)
        (PROJ.4 Params towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.70)
        Accuracy 3m
---
7 Used in Default potsdam region
        (PROJ.4 Params towgs84=606.000,23.000,413.000)
        Default 3-Parameter Transformation

Probably it's not a good idea to automatically
select a datum here (on the other hand most users won't know
which one is the correct datum).

Checking the extent could be done via OGR of course:
ogrinfo -summary strassen-joined.shp strassen-joined
INFO: Open of `strassen-joined.shp'
using driver `ESRI Shapefile' successful.
Layer name: strassen-joined
Geometry: Line String
Feature Count: 12323
Extent: (3427000.290000, 5787594.240000) - (3444004.000000, 5800876.470000)
[...]

(maybe ogrinfo should report in LatLong as well such as gdalinfo does)

The data extent matches:
'3 Used in Germany (West - North - 52d20'N to 55d00'N)'

As usual it's difficult to find out which datum the vendor used
when preparing the data set (hi Intevation!).

>Checking the name, reports:
>grep -i Hauptdreiec ~/grass57/lib/gis/*
>~/grass57/lib/gis/datum.table:potsdam "Deutsches_Hauptdreiecksnetz"
>bessel dx=606.0 dy=23.0 dz=413.0
>
>Since ESRI uses the (gramatically wrong) 'D_Deutsche_Hauptdreiecksnetz',
>GRASS fails to
>find the values. So far no surprise.

But note that the warning message from GRASS used the correct format. I
believe this is due to some good work by Frank in the OSRMorphFromESRI()
from function which is called by GRASS to fix things up before importing
the projection information, so I don't think there is a problem here
really.

Agreed.

>To overcome the problem I would like to add 'D_Deutsche_Hauptdreiecksnetz'
>in lib/gis/datum.table.

The 'correct' way to do this (although it's really a bit of a hack) is to
add two lines to the static const char *papszDatumEquiv in
lib/proj/convert.c and recompile. The first one should be the 'wrong'
value that has occured in your data, and the second one the value in the
GRASS datum.table. I did this for the New Zealand 2000 datum that Hamish
was having some trouble with. Probably the even more correct way to fix it
is to ask Frank to add some code to OGR to fix it.

But as I said I don't think it's necessary here---please check.

Still we have two choices here, 'Whole Country' and 'North-West Germany'
(am I right?). Not sure what to do.
Probably there is no solution for this case (?).

Markus