Coordinate information not stored in tif with native driver

[jipclaassens]

OVSRV6:
C:\ProjDir\Jip\_Github_issues\LandUseModelling_GH648\lus_demo\cfg\demo.dms

/Current_Situation/Current_landuse/Predominant

export data > native driver

native driver tfw content:

100
0
0
-100
50
624950

gdal driver tfw content:

100.0000000000
0.0000000000
0.0000000000
-100.0000000000
100.0000000000
624900.0000000000

Difference means a slight different in output, grids are not aligned.

And opening in ArcGIS shows that the tif from the native driver does not have a spatial reference:

[eoudejans]

Likely a bug on the gdal grid dataset creation side, in gdal_base.cpp line 961:

DPoint offset = ((unit_projection) ? unit_projection->Offset() : DPoint()) + gridBegin * factor + 0.5 * f2;

This means that the offset being used for setting the affine transformation is the middle of the upper-left pixel. This is in contrast to what is expected, namely the upper left coordinate of the upper left pixel:
https://gdal.org/tutorials/geotransforms_tut.html#introduction-to-geotransforms

[eoudejans]

The coordinate information is indeed not stored using the native tif driver at this moment. This can be implemented however.

[eoudejans]

Coordinate information can be stored in GeoTIFF tags:
http://geotiff.maptools.org/spec/geotiff6.html#6.3.3.2

To support the full set of GeoTIFF tags to construct a gdal spatial reference see gdal gtiff driver implementation gt_wkt_srs.cpp:

https://github.com/OSGeo/gdal/blob/master/frmts/gtiff/gt_wkt_srs.cpp#L510
To
https://github.com/OSGeo/gdal/blob/master/frmts/gtiff/gt_wkt_srs.cpp#L1802

[eoudejans]

The gt_wkt_srs.h header is not exposed through gdal installation via vcpkg.

[eoudejans]

Implement a mimimal version of gt_wkt_srs.h using geotiff header: #include "geotiff.h" for writing srs information

[eoudejans]

Mimimal required keys set:

GTIFKeySet(psGTIF, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsArea);
GTIFSetVersionNumbers(psGTIF, GEOTIFF_SPEC_1_1_VERSION, GEOTIFF_SPEC_1_1_KEY_REVISION, GEOTIFF_SPEC_1_1_MINOR_REVISION);
GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1, ModelTypeProjected);
GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
GTIFKeySet(psGTIF, ProjLinearUnitsInterpCorrectGeoKey, TYPE_SHORT, 1, static_cast(1));
GTIFKeySet(psGTIF, GeogLinearUnitsGeoKey, TYPE_SHORT, 1, nUOMLengthCode);
GTIFKeySet(psGTIF, GeogLinearUnitSizeGeoKey, TYPE_DOUBLE, 1, dfLinearUOM);
GTIFKeySet(psGTIF, ProjLinearUnitsGeoKey, TYPE_SHORT, 1, nUOMLengthCode);

[eoudejans]

According to ESRI a TFW file should contain:

Line 1: x-dimension of a pixel in map units
Line 2: rotation parameter
Line 3: rotation parameter
Line 4: NEGATIVE of y-dimension of a pixel in map units
Line 5: x-coordinate of center of upper left pixel
Line 6: y-coordinate of center of upper left pixel

However the SetGeoTranform of gdal is called in gdal_base.cpp to set the geotransform which requires the upper left coordinate for the xy-offset and not the center xy of the upper left pixel as tfw standard specifies:

GT(0) x-coordinate of the upper-left corner of the upper-left pixel.
GT(1) w-e pixel resolution / pixel width.
GT(2) row rotation (typically zero).
GT(3) y-coordinate of the upper-left corner of the upper-left pixel.
GT(4) column rotation (typically zero).
GT(5) n-s pixel resolution / pixel height (negative value for a north-up image).

Adhering to this definition both tfw files are now identical:

tfw native tif:
100
0
0
-100
50
624950

tfw gdal:
100.0000000000
0.0000000000
0.0000000000
-100.0000000000
50.0000000000
624950.0000000000

[eoudejans]

The most urgent part of this issue has been resolved, namely the incorrect tfw file generation. Note that this could have added to translated by +0.5*xres and +0.5-yres. The implementation route of storing coordinate information in the tif header ie complying to the geotiff standard will take more time and will be implemented later.

Note that for implementation the preferred route is to use the PROJ library to provide the necessary projection parameters for storage.