OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.78k stars 2.5k forks source link

gdal_translate RPC adjustment incorrect #10600

Closed dugalh closed 2 weeks ago

dugalh commented 3 weeks ago

What is the bug?

I think gdal_translate does not account for RPC pixel coordinates being relative to the pixel center when adjusting RPC line/sample offsets. The resulting error increases with resolution scaling.

Steps to reproduce the issue

Using this image.

# orthorectify qb2_basic1b.tif
gdalwarp -r bilinear -rpc qb2_basic1b.tif src_rpc_warp.tif

# downsample qb2_basic1b.tif with gdal_translate
gdal_translate -r average -outsize 10% 10% qb2_basic1b.tif gdal_ds.tif

# orthorectify gdal_ds.tif
gdalwarp -r bilinear -rpc gdal_ds.tif gdal_ds_rpc_warp.tif

I see a ~0.5 pixel offset between gdal_ds_rpc_warp.tif and src_rpc_warp.tif if I view them in a GIS.

If I reproduce the downsampled image with Rasterio but using what I think is the correct RPC adjustment:

import rasterio as rio
from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT

# downsample and adjust qb2_basic1b.tif with rasterio
with rio.open('qb2_basic1b.tif', 'r') as src_im:
    scale = 10
    src_win = (0, 0, src_im.width, src_im.height)
    rpcs = src_im.rpcs

    # downsample
    transform = rio.Affine.translation(src_win[0], src_win[1]) * rio.Affine.scale(scale)
    with WarpedVRT(
        src_im,
        transform=transform,
        width=int(src_win[2] / scale),
        height=int(src_win[3] / scale),
        resampling=Resampling.average,
    ) as ds_im:
        array = ds_im.read()

# adjust RPCs (convert center to upper-left convention, scale, then convert back from upper-left
# to center convention)
rpcs.line_off = (rpcs.line_off - src_win[1] + 0.5) / scale - 0.5
rpcs.samp_off = (rpcs.samp_off - src_win[0] + 0.5) / scale - 0.5
rpcs.line_scale /= scale
rpcs.samp_scale /= scale

# write to destination
profile = dict(
    driver='GTiff',
    count=array.shape[0],
    width=array.shape[2],
    height=array.shape[1],
    dtype=array.dtype,
)
with rio.open('user_ds.tif', 'w', **profile, rpcs=rpcs) as dst_im:
    dst_im.write(array)

And orthorectify with gdal_warp:

gdalwarp -r bilinear -rpc user_ds.tif user_ds_rpc_warp.tif

Then user_ds_rpc_warp.tif overlays src_rpc_warp.tif correctly.

Versions and provenance

conda-forge GDAL 3.9.1, released 2024/06/22

Additional context

No response

jratike80 commented 3 weeks ago

Are RPC pixel coordinates always tied to pixel centers, even when the TIFF file has area_or_point=areametadata like qb2_basic1b.tif has?

rouault commented 3 weeks ago

Are RPC pixel coordinates always tied to pixel centers, even when the TIFF file has area_or_point=areametadata like qb2_basic1b.tif has?

yes, at least that's how GDAL implements it. The suggested fix by the OP is correct, and I'm just about to submit a PR with it.