ogr2ogr and gdalwarp -gcp -tps produce different results

[ChrisHSandN]

Expected behavior and actual behavior.

The results of a GCP reprojection seems to be different between the output of gdal_translate & gdalwarp and ogr2ogr

Steps to reproduce the problem.

Given the un-georeferenced TIFF tif_with_cross.zip:

Embed GCP points:

gdal_translate -a_srs 'epsg:4326' \
-gcp '6888' '7412' '-3.8737860398466' '53.21122540405' \
-gcp '2568' '9838' '-3.8477240045557' '53.20254820967' \
-gcp '4808' '9519' '-3.8565346428232' '53.209061881964' \
-gcp '7578' '14429' '-3.8380049721391' '53.227666283207' \
-gcp '3472' '4984' '-3.8728875599878' '53.196744158035' \
-gcp '4751' '5135' '-3.8797033854097' '53.200414741693' \
tif_with_cross.tif tif_with_cross_gdal_translate_gcp.tif

gdapwarp the image:

gdalwarp -tps -t_srs 'epsg:4326' \
tif_with_cross_gdal_translate_gcp.tif \
tif_with_cross_gdal_translate_gcp_gdalwarp.tif

Opening the output of gdalwarp tif_with_cross_gdal_translate_gcp_gdalwarp.tif gives the following result:

The cross on the image is located at pixel 7453,1955, putting this coordinate through ogr2ogr using the same GCP and TPS option:

echo '{"type": "FeatureCollection","features": [{"type": "Feature","geometry": {"coordinates": [7453,1955],"type": "Point"}}]}' | \
ogr2ogr -f GeoJSON \
-gcp '6888' '7412' '-3.8737860398466' '53.21122540405' \
-gcp '2568' '9838' '-3.8477240045557' '53.20254820967' \
-gcp '4808' '9519' '-3.8565346428232' '53.209061881964' \
-gcp '7578' '14429' '-3.8380049721391' '53.227666283207' \
-gcp '3472' '4984' '-3.8728875599878' '53.196744158035' \
-gcp '4751' '5135' '-3.8797033854097' '53.200414741693' \
-tps /vsistdout/ /vsistdin/

Produces:

{
  "type": "FeatureCollection",
  "features": [
	{
	  "type": "Feature",
	  "properties": {},
	  "geometry": {
		"type": "Point",
		"coordinates": [
		  -3.905243393601844,
		  53.20249961030936
		]
	  }
	}
  ]
}

Opening this GeoJSON shows the point is significantly different (~330m) to the location of the cross on the GeoTIFF:

The cross on the GeoTIFF occurs around -3.91015603,53.20182776

Operating system

Rocky Linux 8.6

GDAL version and provenance

GDAL 3.7.2, released 2023/09/05

[jratike80]

I wonder if the y coordinates in the raster case are counted from the top-left corner but in the vector case from the lower-left corner.

[rouault]

The core of the issue is that the TPS (and probably GCP) transformers are not invertible, and for gdalwarp we use the inverse direction to go from a georeferenced location to the source pixel (column, line) coordinate

$ echo 7453 1955 | gdaltransform -tps tif_with_cross_gdal_translate_gcp.tif
-3.90524339360184 53.2024996103094 0

$ echo 7453 1955 | gdaltransform -tps tif_with_cross_gdal_translate_gcp.tif | gdaltransform -i -tps tif_with_cross_gdal_translate_gcp.tif
7223.09557903712 2642.28417742739 0

[ChrisHSandN]

Thanks for your quick reply 👍

I am not sure I quite follow though; both the gdalwarp and ogr2ogr transformations in my example are in the same direction of pixel -> georeferenced.

I am not trying to reverse the transformation, simply expecting, the same transformation done by the 2 utilities to produce the same results for the same input.

To put it another way, if I gave ogr2ogr a set of points for every pixel in the image and then reconstructed them into a raster, should the result not be exactly the same as produced by gdalwarp?

[rouault]

I am not trying to reverse the transformation

That's what you believe ;-) But the way gdalwarp works is that it goes from each target pixel and computes the source pixel coordinate, and that involves the inverse direction of the TPS transformer, which doesn't perfectly roundtrip with the forward way

[ChrisHSandN]

But but, surely one of these results must be inherently “wrong”?!

Given the GCP points and the TPS algorithm, is it not the case that one of these results must be the actual solution and the other is 330m incorrect at that location?

Accepting the inevitable however, what can I do project pixel points in the same way as gdalwarp?
I have known pixel points I wish to overlay on the georeferenced TIFF. Is there some way to reproject this data so that it lines up exactly with the raster produced by gdalwarp?

[rouault]

Working on implementing an exact inverse transformer

[rouault]

With the fix in #8573, we now get much better reversibility

$ echo 7453 1955 | gdaltransform -tps tif_with_cross_gdal_translate_gcp.tif  | gdaltransform -I -tps tif_with_cross_gdal_translate_gcp.tif
7452.99999823652 1955.00000325964 0

[ChrisHSandN]

Perfect, that sounds much more accurate; many thanks for your quick fix. 😃

I will close this issue off now and await the release of 3.8.0.