Grid based transformation with ogr2ogr

classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

Grid based transformation with ogr2ogr

Uwe Schmitz-2
Hi folks,

as a solution for the upcoming datum shift in Germany
(from DHDN to ETRS89[=WGS84]) I'm currently experimenting
with grid based datum shifts.

After fiddling around a little bit, I accomplished to
produce a sample grid in NTv2-Format called DHDN_ETRS89.gsb.

Cs2cs works well with it. The following command-line:

cs2cs +init=epsg:31466 +nadgrids=./DHDN_ETRS89.gsb \
   +to +proj=utm +datum=WGS84 +zone=32

gives me correct transformation results.

Now I'm trying to use ogr2ogr to transform some shape-files.
I'm currently using FWTools1.0.0b2 under WinXP (also
verified with GDAL1.3.1 under Unix) and my command-line
looks like:

ogr2ogr -f "ESRI Shapefile" \
   -s_srs "+init=epsg:31466 +nadgrids=DHDN_ETRS89.gsb" \
   -t_srs "+proj=utm +datum=WGS84 +zone=32" \
   out\GEW01_L.shp in\GEW01_L.shp

If I set PROJ_DEBUG there are only 2 lines of
output:

pj_open_lib(proj_def.dat): call fopen(C:\PROGRA~1\FWTOOL~1.0A6\proj_lib/proj_def
.dat) - succeeded
pj_open_lib(epsg): call fopen(C:\PROGRA~1\FWTOOL~1.0A6\proj_lib/epsg) - succeede
d

Therefore I fear that the grid isn't used at all.
I've also used some other -s_srs/-t_srs combination
with no success

Also I'm wondering why the program doesn't complain, if
I specify a non existant file for +nadgrids.

Thanks for any hints.
Uwe

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Grid based transformation with ogr2ogr

Frank Warmerdam
[hidden email] wrote:

> Hi folks,
>
> as a solution for the upcoming datum shift in Germany
> (from DHDN to ETRS89[=WGS84]) I'm currently experimenting
> with grid based datum shifts.
>
> After fiddling around a little bit, I accomplished to
> produce a sample grid in NTv2-Format called DHDN_ETRS89.gsb.
>
> Cs2cs works well with it. The following command-line:
>
> cs2cs +init=epsg:31466 +nadgrids=./DHDN_ETRS89.gsb \
>    +to +proj=utm +datum=WGS84 +zone=32
>
> gives me correct transformation results.
>
> Now I'm trying to use ogr2ogr to transform some shape-files.
> I'm currently using FWTools1.0.0b2 under WinXP (also
> verified with GDAL1.3.1 under Unix) and my command-line
> looks like:
>
> ogr2ogr -f "ESRI Shapefile" \
>    -s_srs "+init=epsg:31466 +nadgrids=DHDN_ETRS89.gsb" \
>    -t_srs "+proj=utm +datum=WGS84 +zone=32" \
>    out\GEW01_L.shp in\GEW01_L.shp
>
> If I set PROJ_DEBUG there are only 2 lines of
> output:
>
> pj_open_lib(proj_def.dat): call fopen(C:\PROGRA~1\FWTOOL~1.0A6\proj_lib/proj_def
> .dat) - succeeded
> pj_open_lib(epsg): call fopen(C:\PROGRA~1\FWTOOL~1.0A6\proj_lib/epsg) - succeede
> d
>
> Therefore I fear that the grid isn't used at all.
> I've also used some other -s_srs/-t_srs combination
> with no success

Uwe,

Unfortunately gdalwarp uses OGC Well Known Text as the internal coordinate
system representation and that does not include any support for grid files.
NAD27/NAD83 grid shift files are handled as a special case.

So the PROJ.4 representation is parsed and converted to WKT, and then
converted back to PROJ.4 format to do the actual coordinate transformation.
But by then the grid list is lost.

Unfortunately, there is no current way to accomplish what you want with
gdalwarp.  I would encourage you to file a bug report.  I think I may
make an non-standard extension or two to WKT to allow transporting some
of these values I need for PROJ.4.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, [hidden email]
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGF, http://osgeo.org

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

RE: Grid based transformation with ogr2ogr

Martin Daly
In reply to this post by Uwe Schmitz-2
> Unfortunately, there is no current way to accomplish what you
> want with
> gdalwarp.  I would encourage you to file a bug report.  I think I may
> make an non-standard extension or two to WKT to allow
> transporting some
> of these values I need for PROJ.4.

Stop right there, you and your "non-standard extension or two to WKT".

Can you use the "AUTHORITY" WKT construct (<authority> =
AUTHORITY["<name>", "<code>"]), e.g. something like [AUTHORITY="PROJ.4",
"PROJ.4 definition goes here"]?

Regards,
Martin Daly,
Technical Director,
Cadcorp
Computer Aided Development Corporation Ltd.
1, Heathcock Court, LONDON, WC2R 0NT
Tel: +44 (0) 20 7836 0875
Fax: +44 (0) 20 7836 0881
mailto:[hidden email]
http://www.cadcorp.com

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Grid based transformation with ogr2ogr

Frank Warmerdam
Martin Daly wrote:
> Stop right there, you and your "non-standard extension or two to WKT".
 >
> Can you use the "AUTHORITY" WKT construct (<authority> =
> AUTHORITY["<name>", "<code>"]), e.g. something like [AUTHORITY="PROJ.4",
> "PROJ.4 definition goes here"]?

Martin,

Can we have more than one authority entry for a single object?  The logical
place for encoding special information in an authority node would be on
the DATUM but the datum already often has a meaningful relationship to an
EPSG code which I am loath to discard.

Overall, I'm not too keen on embedding a lot of special information in
the authority field.  I'd almost be keener on a special hack to put it into
the name.  I really wish WKT had a "comments" construct to dump this sort
of stuff into.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, [hidden email]
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGF, http://osgeo.org

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

RE: Grid based transformation with ogr2ogr

Martin Daly
In reply to this post by Uwe Schmitz-2
> Can we have more than one authority entry for a single
> object?  

The EBNF only allows one AUTHORITY code.

> The logical
> place for encoding special information in an authority node
> would be on
> the DATUM but the datum already often has a meaningful
> relationship to an
> EPSG code which I am loath to discard.
> Overall, I'm not too keen on embedding a lot of special information in
> the authority field.  I'd almost be keener on a special hack
> to put it into
> the name.  I really wish WKT had a "comments" construct to
> dump this sort
> of stuff into.

Agreed.  What about having both, e.g.
AUTHORITY["EPSG/PROJ.4","1234/PROJ.4 definition"]?  I'd be willing to
handle something like this in our WKT parsing.

M

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

AW: Grid based transformation with ogr2ogr

Uwe Schmitz-2
In reply to this post by Uwe Schmitz-2
Hi Frank,

>
> Uwe,
>
> Unfortunately gdalwarp uses OGC Well Known Text as the
> internal coordinate
> system representation and that does not include any support
> for grid files.
> NAD27/NAD83 grid shift files are handled as a special case.
>
> So the PROJ.4 representation is parsed and converted to WKT, and then
> converted back to PROJ.4 format to do the actual coordinate
> transformation.
> But by then the grid list is lost.
>
> Unfortunately, there is no current way to accomplish what you
> want with
> gdalwarp.  I would encourage you to file a bug report.  I think I may
> make an non-standard extension or two to WKT to allow
> transporting some
> of these values I need for PROJ.4.
>

ok, I see the problem. I will file a bug report
as soon as I manage to log in (It seems that I've
forgotten my password... and I'm waiting for the
mail to set a new one).

Meanwhile... is it possible to do the datum-shift
with the +towgs84=x,y,z switch, or is it lost in
the same way? And if so, where has +towgs84 to
be placed? In the -s_srs or the -t_srs switch?

Many thanks in advance
Uwe

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: AW: Grid based transformation with ogr2ogr

Frank Warmerdam
[hidden email] wrote:
> ok, I see the problem. I will file a bug report
> as soon as I manage to log in (It seems that I've
> forgotten my password... and I'm waiting for the
> mail to set a new one).
>
> Meanwhile... is it possible to do the datum-shift
> with the +towgs84=x,y,z switch, or is it lost in
> the same way? And if so, where has +towgs84 to
> be placed? In the -s_srs or the -t_srs switch?

Uwe,

The +towgs84 parameter will be preserved as a TOWGS84[] entry
for the datum in WKT.  So that works fine.  In your case it would
belong in the -s_srs since it relates the DHDN to WGS84.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, [hidden email]
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGF, http://osgeo.org

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Grid based transformation with ogr2ogr

Uwe Schmitz-2
In reply to this post by Uwe Schmitz-2
Frank,

> Uwe,
>
> The +towgs84 parameter will be preserved as a TOWGS84[] entry
> for the datum in WKT.  So that works fine.  In your case it would
> belong in the -s_srs since it relates the DHDN to WGS84.
>

forgive me if I get on your nerves but...

I tried the following (this time I use GML for better reading,
with shape files I gain the same result):

ogr2ogr -f "GML" -s_srs "+init=epsg:31467 +towgs84=631,23,451" \
   -t_srs "+proj=utm +datum=WGS84 +zone=32" trf.xml org.xml

and for:
...<gml:coordinates>3396682.3999999999,5626850.4000000004</gml:coordinates>...

in org.xml I get:
...<gml:coordinates>396651.68538503139,5625036.6721786009</gml:coordinates>...

in trf.xml.

So far so good but:
cs2cs +init=epsg:31467 +towgs84=631,23,451 \
   +to +proj=utm +datum=WGS84 +zone=32<CR>
3396682.3999999999 5626850.4000000004 0.0<CR>

gives me:
396648.50       5625041.57 46.88

IMO cs2cs is correct.
What's going wrong?

Best wishes
Uwe

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Grid based transformation with ogr2ogr

Frank Warmerdam
[hidden email] wrote:
> forgive me if I get on your nerves but...
>
> I tried the following (this time I use GML for better reading,
> with shape files I gain the same result):
>
> ogr2ogr -f "GML" -s_srs "+init=epsg:31467 +towgs84=631,23,451" \
>    -t_srs "+proj=utm +datum=WGS84 +zone=32" trf.xml org.xml

> and for:
> ...<gml:coordinates>3396682.3999999999,5626850.4000000004</gml:coordinates>...
>
> in org.xml I get:
> ...<gml:coordinates>396651.68538503139,5625036.6721786009</gml:coordinates>...

Uwe,

I must appologising for not noticing/commenting on another issue.  If
you expand +init=epsg:31467 it looks like:

   +proj=tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0
      +ellps=bessel +datum=potsdam +units=m +no_defs

The +no_defs part actually means no defaults, but it is implemented
by ignoring all items that come after it.  So I believe that the
+towgs84= clause is being ignored by one or both of ogr2ogr and
cs2cs.  To ensure you control the settings you should provide the
definition directly instead of using the +init mechanism.

eg.
    +proj=tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0
      +ellps=bessel +towgs84=631,23,451 +units=m +no_defs

I also discovered that you do *not* want to use +datum=potsdam or
else the internal PROJ.4 towgs84 definition for potsdam will take over.

I must confess this stuff isn't as obvious as it ought to be!


> in trf.xml.
>
> So far so good but:
> cs2cs +init=epsg:31467 +towgs84=631,23,451 \
>    +to +proj=utm +datum=WGS84 +zone=32<CR>
> 3396682.3999999999 5626850.4000000004 0.0<CR>
>
> gives me:
> 396648.50       5625041.57 46.88
>
> IMO cs2cs is correct.
> What's going wrong?

Right.  I have confirmed though that:

testepsg -t '+proj=tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0
+ellps=bessel +units=m +towgs84=631.0,23.0,451.0' '+proj=utm +datum=WGS84
+zone=32' 3396682.3999999999 5626850.4000000004

gives the same result.  This uses the OGRSpatialReference layer, similar
to what ogr2ogr will do.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, [hidden email]
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGF, http://osgeo.org

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Grid based transformation with ogr2ogr

Uwe Schmitz-2
In reply to this post by Uwe Schmitz-2
Frank,

> Uwe,
>
> I must appologising for not noticing/commenting on another issue.  If
> you expand +init=epsg:31467 it looks like:
>
>    +proj=tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0
>       +ellps=bessel +datum=potsdam +units=m +no_defs
>
> The +no_defs part actually means no defaults, but it is implemented
> by ignoring all items that come after it.  So I believe that the
> +towgs84= clause is being ignored by one or both of ogr2ogr and
> cs2cs.  To ensure you control the settings you should provide the
> definition directly instead of using the +init mechanism.
>
> eg.
>     +proj=tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0
>       +ellps=bessel +towgs84=631,23,451 +units=m +no_defs
>
> I also discovered that you do *not* want to use +datum=potsdam or
> else the internal PROJ.4 towgs84 definition for potsdam will
> take over.
>
> I must confess this stuff isn't as obvious as it ought to be!
>

Ok, now it works as expected! Our discussion clarifies
some things for me.

Many thanks for your help!
Uwe

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Grid based transformation with ogr2ogr

mhw-at-yg
In reply to this post by Frank Warmerdam
Hi Uwe,Frank,

>> as a solution for the upcoming datum shift in Germany (from DHDN to
>>  ETRS89[=WGS84]) I'm currently experimenting with grid based datum
>> shifts.
...

> Unfortunately gdalwarp uses OGC Well Known Text as the internal
> coordinate system representation and that does not include any
> support for grid files. NAD27/NAD83 grid shift files are handled as a
> special case.
>
> So the PROJ.4 representation is parsed and converted to WKT, and then
>  converted back to PROJ.4 format to do the actual coordinate
> transformation. But by then the grid list is lost.
>
> Unfortunately, there is no current way to accomplish what you want
> with gdalwarp.  I would encourage you to file a bug report.  I think
> I may make an non-standard extension or two to WKT to allow
> transporting some of these values I need for PROJ.4.


To my eye this looks related to my earlier proposal to develop an ogr
module to apply the Canadian Correction Matrices (CORMAT) to vector
data. Is that correct? If these are similar problems perhaps there is an
opportunity to share resources.

cheers,

--
matt wilkie
--------------------------------------------
Geographic Information,
Information Management and Technology,
Yukon Department of Environment
10 Burns Road * Whitehorse, Yukon * Y1A 4Y9
867-667-8133 Tel * 867-393-7003 Fax
http://environmentyukon.gov.yk.ca/geomatics/
--------------------------------------------

_______________________________________________
Gdal-dev mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/gdal-dev
-Matt