gdal_retile moves tiles a few KM from where their source GeoTIFF's footprint was

[marklit]

What is the bug?

After running gdal_retile, the resulting smaller tiles are a few KM away from their parent's footprint. Here the dark tiles should be over the airport but they're offset by a few KM. The tiles produce a smaller footprint overall as well. The source image covered the entire Airport's footprint. They're also rotated incorrectly.

Steps to reproduce the issue

$ wget https://umbra-open-data-catalog.s3.amazonaws.com/sar-data/tasks/ad%20hoc/Haneda_Japan/0028c361-193c-49b7-ad59-808033d5dc74/2024-02-27-01-18-51_UMBRA-07/2024-02-27-01-18-51_UMBRA-07_CSI.tif
$ mkdir -p tiles
$ gdal_retile.py \
    -s_srs "EPSG:4326" \
    -ps 4096 4096 \
    -targetDir tiles/ \
    2024-02-27-01-18-51_UMBRA-07_CSI.tif

Versions and provenance

GDAL 3.4.1, released 2021/12/27

I'm using this version as it's the latest Ubuntu for Windows Package Manager supports.

Additional context

➫ gdalinfo 2024-02-27-01-18-51_UMBRA-07_CSI.tif

Driver: GTiff/GeoTIFF
Files: 2024-02-27-01-18-51_UMBRA-07_CSI.tif
Size is 16001, 16001
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
GeoTransform =
  139.8071184948572, 9.877606409126224e-07, -3.301795633394074e-06
  35.576945271941, -2.698390246796095e-06, -8.072470111920523e-07
Metadata:
  AREA_OR_POINT=Area
  COLLECT_ID=898410b9-67ab-4f7d-b6d7-688f53fd3d38
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
  LAYOUT=COG
RPC Metadata:
  ERR_BIAS=-1
  ERR_RAND=-1
  HEIGHT_OFF=36.5424368376879
  HEIGHT_SCALE=1414.21349961605
  LAT_OFF=35.5489016148089
  LAT_SCALE=0.0138154544715561
  LINE_DEN_COEFF=1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  LINE_NUM_COEFF=4.44272318597276e-07 -0.557824708958129 -0.166878164344482 -0.812869251845968 9.58393511233691e-05 -0.00405020664411899 0.000720521045028307 9.10354435475097e-05 0.000106030390123858 -0.00206906795942957 7.7329377257039e-06 -7.11815349487315e-07 -2.31456828350157e-07 -4.16516985215374e-05 -7.6414197251492e-07 -6.60516165613977e-07 7.46360408355076e-06 -2.36426705765341e-05 -5.45919285636177e-06 -3.88591220272253e-05
  LINE_OFF=8000
  LINE_SCALE=8424.3261229398
  LONG_OFF=139.788605057899
  LONG_SCALE=0.0169048199998421
  SAMP_DEN_COEFF=1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  SAMP_NUM_COEFF=-2.96052985813837e-07 0.231574776742832 -0.774086452562128 0.589086057957499 -3.98269216991867e-05 0.00281093446760478 -0.000704678120832472 -0.000142748949521291 -7.78699184718959e-05 0.00145412159792637 -5.42592529219741e-06 3.03039659816264e-07 9.57745474447393e-07 2.92716769311354e-05 -1.22154126364717e-06 -2.06105745822261e-06 -6.51652196885433e-06 1.81791221177082e-05 5.41548640560391e-06 3.4055873360515e-05
  SAMP_OFF=8000
  SAMP_SCALE=6070.75035808271
Corner Coordinates:
Upper Left  ( 139.8071185,  35.5769453) (139d48'25.63"E, 35d34'37.00"N)
Lower Left  ( 139.7542865,  35.5640285) (139d45'15.43"E, 35d33'50.50"N)
Upper Right ( 139.8229237,  35.5337683) (139d49'22.53"E, 35d32' 1.57"N)
Lower Right ( 139.7700916,  35.5208516) (139d46'12.33"E, 35d31'15.07"N)
Center      ( 139.7886051,  35.5488984) (139d47'18.98"E, 35d32'56.03"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Undefined
  NoData Value=0
Band 2 Block=512x512 Type=Byte, ColorInterp=Undefined
  NoData Value=0
Band 3 Block=512x512 Type=Byte, ColorInterp=Undefined
  NoData Value=0

[rouault]

gdal_retile is neither ready for source files that have rotation terms in their geotransform, nor RPC. Currently you would probably be better implementing some sort of gdal_retile logic by scripting gdal_translate, which should be more ready for the task

[marklit]

Understood. I'll look into that.

What do you mean by RPC?

[rouault]

What do you mean by RPC?

The ugly set of stuff below the "RPC Metadata:" line in your gdalinfo output ;-) cf https://en.wikipedia.org/wiki/Rational_polynomial_coefficient and https://gdal.org/development/rfc/rfc22_rpc.html

[marklit]

I asked Umbra for advice on a set of gdal parameters https://twitter.com/marklit82/status/1809151161825427959

[marklit]

A single call to gdalwarp fixed it.

$ wget https://umbra-open-data-catalog.s3.amazonaws.com/sar-data/tasks/ad%20hoc/Haneda_Japan/0028c361-193c-49b7-ad59-808033d5dc74/2024-02-27-01-18-51_UMBRA-07/2024-02-27-01-18-51_UMBRA-07_CSI.tif

$ gdalwarp -t_srs "EPSG:4326" 2024-02-27-01-18-51_UMBRA-07_CSI.tif warped.tif

$ gdal_retile.py \
    -s_srs "EPSG:4326" \
    -ps 4096 4096 \
    -targetDir ./ \
    warped.tif