[gdal-dev] how to use the -ct option in ogr2ogr

[Markus Metz]

On Mon, Feb 8, 2021 at 10:23 AM Karsten Tebling <tebling at masuch.de> wrote:
>
> gdalinfo for NTv2LSBB_LSA.gsb gave the following information:
>
> Size is 256, 264
> Coordinate System is: EPSG 4326
> Origin: 10.5, 53.1
> SYSTEM_F: S4283
> SYSTEM_T: ETRS89
> Upper Left 10.5, 53.1
> Lower Left 10.5, 50.9
> Upper Right 13.3, 53.1
> Lower Right 13.3, 50.9
>
> ogrinfo for punkte_31468.dxf gave the following information:
>
> Geometry: Unknown (any)
> Feature Count: 36
> Extent: (4523806.5, 5758974.7) - (4524134.7, 5759234.0)
> Layer SRS WKT: (unknown)
>
> according to QGIS punkte_31468.dxf should be at 11.6 52.1 which is inside
the extent of the NTv2LSBB_LSA.gsb. after checking the documentation for
the NTv2LSBB_LSA.gsb i noticed i couldn't use it for transformations from
31468 to 25833, that's why i tried it with BETA2007.gsb instead, but it
didn't work either - it gave the same errors.
>
> gdalinfo for BETA2007.gsb gave the following information:
>
> Size is 62, 84
> Coordinate System is: EPSG 4326
> Origin: 5.4, 55.4
> SYSTEM_F: DHDN90
> SYSTEM_T: ETRS89
> Upper Left 5.4, 55.4
> Lower Left 5.4, 47.0
> Upper Right 15.8, 55.4
> Lower Right 15.6, 47.0
>
> my adjusted code for BETA2007.gsb looks like:
>
> ogr2ogr -s_srs "EPSG:31468" -t_srs "EPSG:25833" -dim XYZ -ct
"+proj=pipeline +step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000
+y_0=0 +ellps=bessel +step +proj=hgridshift +grids=BETA2007.gsb +step
+proj=utm +zone=33 +ellps=GRS80" -f DXF punkte_25833.dxf punkte_31468.dxf
>
> the points i used to test this were digitized around Magdeburg,
Saxony-Anhalt in EPSG:31468. according to aerial data and openstreetmap,
the position and coordinate system are correct when i check it in QGIS and
switch to different coordinate systems.

Wit PROJ 7.2.1 and
projinfo -o PROJ -s "EPSG:31468" -t "EPSG:25833"
I get
Operation No. 1:

unknown id, Inverse of 3-degree Gauss-Kruger zone 4 + DHDN to ETRS89 (8) +
UTM zone 33N, 0.9 m, Germany - onshore - states of Baden-Wurtemberg,
Bayern, Berlin, Brandenburg, Bremen, Hamburg, Hessen,
Mecklenburg-Vorpommern, Niedersachsen, Nordrhein-Westfalen,
Rheinland-Pfalz, Saarland, Sachsen, Sachsen-Anhalt, Schleswig-Holstein,
Thuringen.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
        +ellps=bessel
  +step +proj=hgridshift +grids=de_adv_BETA2007.tif
  +step +proj=utm +zone=33 +ellps=GRS80

-------------------------------------
Operation No. 2:

unknown id, Inverse of 3-degree Gauss-Kruger zone 4 + DHDN to ETRS89 (2) +
UTM zone 33N, 3 m, Germany - states of former West Germany onshore -
Baden-Wurtemberg, Bayern, Bremen, Hamburg, Hessen, Niedersachsen,
Nordrhein-Westfalen, Rheinland-Pfalz, Saarland, Schleswig-Holstein.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
        +ellps=bessel
  +step +proj=push +v_3
  +step +proj=cart +ellps=bessel
  +step +proj=helmert +x=598.1 +y=73.7 +z=418.2 +rx=0.202 +ry=0.045
+rz=-2.455
        +s=6.7 +convention=position_vector
  +step +inv +proj=cart +ellps=GRS80
  +step +proj=pop +v_3
  +step +proj=utm +zone=33 +ellps=GRS80


It seems your custom definitions of the coordinate operation are missing
the axis swap step "+step +proj=axisswap +order=2,1".
According to "projinfo -o WKT2:2019 "EPSG:31468", northing is the first
axis in EPSG:31468, thus the need for the axis swap.

If you trust the QGIS coordinate operation, then why do you provide a
custom coordinate operation to ogr2ogr with the -ct option?

>
> Am 06.02.2021 um 22:24 schrieb Markus Metz:
>
>
>
> On Fri, Feb 5, 2021 at 3:40 PM Karsten Tebling <tebling at masuch.de> wrote:
> >
> > Hello,
> >
> > i'm struggling with how to use the -ct option in ogr2ogr and could need
> > some help.
> >
> > i copied a NTv2.gsb to osgeo4w64/share/proj and edited the proj.db by
> > inserting a new row in the grid_transformation-table. using this
> > transformation in qgis works. however when i try to use it in ogr2ogr
> > with the same data i always get "failed to reproject feature (geometry
> > probably out of source or destination srs). no known way to write
> > feature with geometry 'none'". here is my code:
> >
> > ogr2ogr -s_srs "EPSG:31468" -t_srs "EPSG:25833" -dim XYZ -ct
> > "+proj=pipeline +step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1
> > +x_0=4500000 +y_0=0 +ellps=bessel +step +proj=hgridshift
> > +grids=NTv2LSBB_LSA.gsb +step +proj=utm +zone=33 +ellps=GRS80" -f DXF
> > punkte_25833.dxf punkte_31468.dxf
> >
>
> Is the grid NTv2LSBB_LSA.gsb used for hgridshift really covering the full
extent of punkte_31468.dxf ?
> You can check with "gdalinfo NTv2LSBB_LSA.gsb" and "ogrinfo -so
punkte_31468.dxf".
>
> Markus M
>
> > i'm using the latest qgis-ltr and osgeo4w shell
> >
> > greetings
> > _______________________________________________
> > gdal-dev mailing list
> > gdal-dev at lists.osgeo.org
> > https://lists.osgeo.org/mailman/listinfo/gdal-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20210208/567135c8/attachment.html>