gdalwarp with geolocation arrays: Systematic Pixel shift

[piyushrpt]

Hi,

We have been testing out the use of geolocation arrays to geocode swath data and have noticed a systematic one-pixel shift either in X or Y directions based on the orientation of the swath. To verify this I created the following unit test:

1. Create a rotated lat.tif, lon.tif grid. And use either lat.tif or lon.tif as the data to be reprojected itself.
2. Compare the reprojected data with what is expected from the geotransform for the reprojected data.
3. I'm attaching the python script that I used to run the tests.
4. I'm also attaching two text files with the outputs of the python script - one where I let gdalwarp pick output spacing and other where I use an oversampling factor of 4.

In either case, the impact of orientation is clear. There is no shift (within noise) observed when the swath is almost approximately oriented N-S, W-E to begin with.

This could either be an issue in the interpolation or when the geotransform is setup for the reprojected output.

Piyush

[piyushrpt]

Stats when gdalwarp is allowed to pick output spacing

[piyushrpt]

Output when an oversampling factor of 4 is used for output grid

[piyushrpt]

Python test script

[piyushrpt]

Test script with EW, NS aligned grids

[piyushrpt]

Output for aligned grids with oversampling factor of 4

[piyushrpt]

I further simplified the unit test and used grids oriented in EW, NS directions. The only orientation of the geoloc arrays that doesn't result in a systematic shift is WE-NS as expected. If any of the axes is reversed (i.e, orientation of the geoloc arrays) we see a systematic shift of 1 pixel. I'm attaching the simpler test script with output at oversampling factor 4. You will notice the pixel shift of 3.x pixels in Means.

Piyush

[Even Rouault]

Piyush, I doubt I'll be able to look at that on my free time / voluntary basis. From the quality of the bug report and your previous contributions, I believe you might perhaps have the capability to debug this by yourself. If so, you can concentrate on alg/gdalgeoloc.cpp . I believe that with gdalwarp we are taking the reverse transformer (to be confirmed), so the "geox/geoy to pixel/line using backmap." case at line 832. But the root cause is probably in the way the backmap is built in GeoLocGenerateBackMap() (line 213). Perhaps a quick way to decrease the error would be to try to decrease the pixel size of the backmap (line 259)

[piyushrpt]

Hi Even,

Thanks for pointing me to the right sections of the code. Was unsure about what to look at. Will go through this in detail and get back to you.

Piyush

[piyushrpt]

I'm almost convinced through testing that the issue is due to the fact that the mapping indices are truncated using "static_cast<int>" when the appropriate thing to do would be to round the floating point number and then cast it to int. Is that a valid solution to pursue? Do you foresee other unexpected consequences?

[Even Rouault]

Well if rounding instead of truncating improves things in some cases and does not degrade in cases where it used to work, it looks like it might be a valid fix. It might be convenient for you to submit as a pull request against the github clone so as to see if that breaks existing tests (which might then probably just require adjusting their expected results)

[piyushrpt]

Hi Even,

I spent quite some time verifying the backmap generation in both Python and C++. These seem consistent for the aligned orientations that I used for the unit tests above - WENS, EWNS, WESN, EWSN. The outputs of GDALGeoLocTransform seem mostly consistent with input data in all those orientations. Will keep digging.

Piyush

[piyushrpt]

Test script for geolocation arrays - different orientations

[piyushrpt]

Discrepancy plot for WENS. Error is same for all orientations. Left is discrepancy in "y" and Right is discrepancy in "x"

[Even Rouault]

Has been fixed per ​https://github.com/OSGeo/gdal/pull/244 r40391

[warmerdam]

I am seeing significant problems still with inverse transformation with the geolocation transformer. I'm reopening this ticket, and I'll try to prepare a fairly minimalist demonstration of the issue and potentially a fix though the fix is not yet obvious to me.

Note, my usecase is geolocation arrays that represent "fine registration" results for which the backmap approach may not be all that appropriate. But initially I'm just doing simple grids to demonstrate the problem.

[warmerdam]

Files to reproduce current errors

[warmerdam]

to reproduce the problem I'm seeing, unpack "repro2.zip" (attached). This has a reduced resolution mesh described by:

Geolocation:
  LINE_OFFSET=0.0
  LINE_STEP=60.0
  PIXEL_OFFSET=0.0
  PIXEL_STEP=64.0
  SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
  X_BAND=1
  X_DATASET=mesh.tif
  Y_BAND=2
  Y_DATASET=mesh.tif

The mesh.tif just values mapping pixel/line locations to themselves in pixel/line space (though I label it as WGS84).

warmerdam@gdal-tp:~/pl/planet_common/wrk5/repro2$ ./mkadj.py 
[[[    0.    64.   128. ...,  2432.  2496.  2560.]
  [    0.    64.   128. ...,  2432.  2496.  2560.]
  [    0.    64.   128. ...,  2432.  2496.  2560.]
  ..., 
  [    0.    64.   128. ...,  2432.  2496.  2560.]
  [    0.    64.   128. ...,  2432.  2496.  2560.]
  [    0.    64.   128. ...,  2432.  2496.  2560.]]

 [[    0.     0.     0. ...,     0.     0.     0.]
  [   60.    60.    60. ...,    60.    60.    60.]
  [  120.   120.   120. ...,   120.   120.   120.]
  ..., 
  [  960.   960.   960. ...,   960.   960.   960.]
  [ 1020.  1020.  1020. ...,  1020.  1020.  1020.]
  [ 1080.  1080.  1080. ...,  1080.  1080.  1080.]]]

Then, if I use gdaltransform to transform to "geolocation" space I get reasonable values:

vagrant@vagrant:/vagrant/wrk5/repro2$ gdaltransform src.vrt -t_srs WGS84
0 0
0 0 0
64 60 
64 60 0
32 30
32 30 0

But if I do the invert, I get strange results:

vagrant@vagrant:/vagrant/wrk5/repro2$ gdaltransform src.vrt -t_srs WGS84 -i
0 0
32 66.7246475219727 0
64 60
106.216476396658 90.4450774047793 0
32 30
71.1957064920051 66.7246428494803 0

I am vaguely suspecting the problem relates to different step size in X and Y though I have not yet experimented with that, but things are clearly broken in this trivial mapping of values to themselves.

[warmerdam]

an even more trivial case in repro3.zip which is a 10x10 src.vrt, and a 10x10 mesh that maps each pixel to itself. To step size, no asymmetry, no odd numbers.

vagrant@vagrant:/vagrant/wrk5/repro3$ gdaltransform -geoloc src.vrt -t_srs WGS84 -i
0 0 
0.5 0.683993399143219 0
5 5
5.62062359025868 5.43875897040691 0
10 10 
transformation failed.
9 9 
9.31600666046143 9.5 0

The 0.5 error looks to me like it is the FSHIFT/ISHIFT introduced in @piyushrpt's change, but there is still additional error coming from somewhere. I will dig deeper.

[warmerdam]

simplified 10x10 example of the problem

[warmerdam]

I'm coming to the conclusion that:

• the forward transformation is not actually bilinear (influence of the x+1,y+1 corner of the grid cell is ignored as it is approached).
• the backward map is fundamentally broken and a bad idea except perhaps for an initial guess.

I'm contemplating making a new implementation of the geolocation array transformer that is truely bilinear in the forward direction and uses a generic iterative inverse approach (possibly with something like the backward map for initial guesses). I would like the geolocation transformer to be a useful general purpose mesh transformer.

[warmerdam]

correction, the forward transform is bilinear and is working fine. I had mis-parsed the code.

[warmerdam]

FYI, I've backed off work on a new implementation though it is possible I might return to it at some point. Currently this is available for others to work on.

[liporace]

Replying to warmerdam:

an even more trivial case in repro3.zip which is a 10x10 src.vrt, and a 10x10 mesh that maps each pixel to itself. To step size, no asymmetry, no odd numbers.

vagrant@vagrant:/vagrant/wrk5/repro3$ gdaltransform -geoloc src.vrt -t_srs WGS84 -i
0 0 
0.5 0.683993399143219 0
5 5
5.62062359025868 5.43875897040691 0
10 10 
transformation failed.
9 9 
9.31600666046143 9.5 0

The 0.5 error looks to me like it is the FSHIFT/ISHIFT introduced in @piyushrpt's change, but there is still additional error coming from somewhere. I will dig deeper.

I've checked this example against the current (master) implementation, with FSHIFT and ISHIFT set to 0 and an increased OVERSAMPLE_FACTOR and the results seem to be OK:

(gdal) gdaltransform -geoloc src.vrt -t_srs WGS84 -i
Enter column line values separated by space, and press Return.
0 0
0 0 0
5 5
4.99999975127938 4.99999975127938 0
10 10
transformation failed.
9 9
9 9 0

If I decrease OVERSAMPLE_FACTOR then I start to get some differences:

(gdal) gdaltransform -geoloc src.vrt -t_srs WGS84 -i
Enter column line values separated by space, and press Return.
0 0
0.0499999970197678 0.0999999940395355 0
5 5
4.8611108623905 4.86111089918348 0
10 10
transformation failed.
9 9
8.89999961853027 8.95000028610229 0

Maybe this is a matter of setting an appropriate factor?