Local coordinate system to EPSG via affine transformation

Question

I have a raster file with a local coordinate system. It is based on EPSG:31466 but with a specific rotation and offset.

gdalinfo output is this:

I need to transform this gtiff back to 31466, but how?

With FME I can use the Affiner transformer to do this.

Affiner parameter:

gdalinfo of output geotiff from FME:

But I need to know the workflow with gdal.

I'm using gdal 2.3.0.

Answer 1

Since FME writes a Geotransform table into the metadata, you can just take that for your mosaik file:

Transform your file to VRT

gdal_translate -of VRT mosaik.tif mosaik.vrt

Open mosaik.vrt in a text editor and insert a Geotransform line:

<VRTDataset rasterXSize="15000" rasterYSize="25000">
<GeoTransform>  2575426.714515543,  0.011400621271025,  0.004378508269806,  5729751.965613389,  0.004376888527729, -0.011404849098262</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1">
    <ColorInterp>Red</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">mosaik.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="15000" RasterYSize="25000" DataType="Byte" BlockXSize="15000" BlockYSize="64" />
      <SrcRect xOff="0" yOff="0" xSize="15000" ySize="25000" />
      <DstRect xOff="0" yOff="0" xSize="15000" ySize="25000" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

Alternatively, you can create a custom omerc projection around the image center, as described in How to create proj4 definition of local coordinatesystem

+proj=omerc +lat_0=51.6959777875 +lonc=7.0923165808 +alpha=-20.145 +gamma=0 +k=1 +x_0=50692.579 +y_0=81723.458 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs