Gdalwarp orthorectification WorldView-3 not using RPC projection properly

Question

I have a high-resolution WorldView-3 image that I am trying to orthorectify using

gdalwarp -rpc -to "RPC_DEM=dem.tif" input.tif output.tif

However the output is not aligned with Bing or Google maps at all,

I have tried various DEM sources both in WSG84 as GeoID format, but both give bad results and also tried different target EPSGs.

When I manually use the RPC model to orthorectify small tiles of the image by projecting the corners using code from https://github.com/gfacciol/IS18 the tiles do match Bing and Google maps very well. However this is a slow process and results in edges of the tiles not matching exactly. But it shows that the RPC model does contain the correct information to orthorectify the image.

This is the code that I use to orthorectify a tile within the image that does seem to work (even though it seems like a stupid solution):

from pyproj import Proj, transform
from rasterio.windows import Window
import utils # https://github.com/gfacciol/IS18/blob/master/utils.py
outProj = Proj(init='epsg:32629')
inProj = Proj(init='epsg:4326')
x1, y1 = transform(inProj,outProj,aoi['coordinates'][0][0][0], aoi['coordinates'][0][0][1])
z = srtm4.srtm4(lon_center, lat_center) # determines ellipsoidal height at AOI
myrpc = utils.rpc_from_geotiff(image)
x1, y1 = myrpc.projection(lon_upperleft, lat_upperleft, z)
x2, y2 = myrpc.projection(lon_bottomright, lat_bottomright, z)
easting, northing = transform(inProj,outProj, lon_upperleft, lat_upperleft)
easting2, northing2 = transform(outProj,inProj,lon_bottomright, lat_bottomright)
w, h = x2 - x1, y2 - y1
window = Window(x1, y1, w, h)
metadata = ds.meta
metadata.update(transform=rasterio.transform.from_origin(easting, northing, 0.3, 0.3))
metadata.update(width=w, height=h)

with rasterio.open('test_ortho.tif', 'w', **metadata) as dst:
    dst.write_band(1, ds.read(window=window)[0,:,:])

This is the output I get with --config CPL_DEBUG ON:

gdalwarp --config CPL_DEBUG ON  -rpc -to RPC_DEM="/home/ubuntu/SRTM/srtm_WSG_ellipsoid.tif" input.tif output.tif
GDAL: GDALOpen(input.tif, this=0x56001a5b0880) succeeds as GTiff.
MDReaderDigitalGlobe: RPB Filename: input.RPB
GDAL: GDALOpen(/home/ubuntu/SRTM/srtm_WSG_ellipsoid.tif, this=0x56001a5b69e0) succeeds as GTiff.
OGRCT: PROJ >= 4.8.0 features enabled
OGRCT: Using locale-safe proj version
OGRCT: Source: +proj=longlat +datum=WGS84 +no_defs
OGRCT: Target: +proj=longlat +datum=WGS84 +no_defs
GDAL: GDAL_CACHEMAX = 3069 MB
Creating output file that is 18490P x 14370L.
GDAL: GDALDriver::Create(GTiff,output.tif,18490,14370,1,UInt16,0)
Processing input file input.tif.
WARP: Copying metadata from first source to destination dataset
GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=0,0,4061x4083 Dst=0,0,4622x3592
0..GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=4050,0,4112x4030 Dst=4622,0,4623x3592
.10GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=8153,0,2066x3996 Dst=9245,0,2311x3592
..GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=10204,0,2108x4122 Dst=11556,0,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=12263,0,2065x4113 Dst=13867,0,2311x3592
20GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=14322,0,2062x4006 Dst=16178,0,2312x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=0,4025,4052x4117 Dst=0,3592,4622x3593
..30GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=4030,3996,2069x4121 Dst=4622,3592,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=6085,3993,2073x4106 Dst=6933,3592,2312x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=8143,3969,2062x4136 Dst=9245,3592,2311x3593
.40GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=10201,3968,2126x4157 Dst=11556,3592,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=12269,3978,2109x4274 Dst=13867,3592,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=14322,3978,2062x4280 Dst=16178,3592,2312x3593
.50GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=0,8125,1981x4248 Dst=0,7185,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=1966,8105,2065x4124 Dst=2311,7185,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=4024,8097,2065x4138 Dst=4622,7185,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=6084,8093,2061x4148 Dst=6933,7185,2312x3592
60GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=8141,8091,2062x4140 Dst=9245,7185,2311x3592
..GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=10199,8095,2078x4159 Dst=11556,7185,2311x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=12267,8124,2138x4208 Dst=13867,7185,2311x3592
70GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=14344,8229,2040x4158 Dst=16178,7185,2312x3592
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=0,12224,2054x4160 Dst=0,10777,2311x3593
..GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=1967,12223,2143x4161 Dst=2311,10777,2311x3593
80GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=4023,12224,2140x4160 Dst=4622,10777,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=6083,12226,2104x4158 Dst=6933,10777,2312x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=8142,12225,2108x4159 Dst=9245,10777,2311x3593
.90GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=10200,12227,2133x4157 Dst=11556,10777,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=12264,12253,2113x4131 Dst=13867,10777,2311x3593
.GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyShort() Src=14324,12303,2060x4081 Dst=16178,10777,2312x3593
.100 - done.
GDAL: GDALClose(/home/ubuntu/SRTM/srtm_WSG_ellipsoid.tif, this=0x56001a5b69e0)
GDAL: GDALClose(output.tif, this=0x56001a5fa700)
GDAL: GDALClose(input.tif, this=0x56001a5b0880)

This is in GDAL 2.2.3.

In GDAL 3.1.0dev I get an additional debug output (the rest is basically the same):

RPC: Short-circuiting coordinate transformation from DEM SRS to WGS 84 due to apparent nop
GDAL: GDAL_CACHEMAX = 3069 MB
RPC: Using GDALRPCTransformWholeLineWithDEM

Does anyone have an idea what is happening here?

Answer 1

I finally managed to correctly orthorectify my image using Gdal 2.2.2.

I am not 100% sure which made the difference, but there are two things that I did.

Firstly I made sure that the DEM used has the correct datum (WSG84 ellipsoid and not a geoid datum) by applying:

gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=egm08_25.gtx" -t_srs "+proj=longlat +datum=WGS84 +no_def" DEM.tif DEM_WSG84_ellipsoid.tif

In this case I used egm08_25.gtx since my DEM was in EGM2008 geoid height, but for NASA's SRTM it should be EGM96.

Then when warping the image I make sure to add the -et 0 option as suggested here: https://lists.osgeo.org/pipermail/gdal-dev/2016-July/044907.html which might be the reason it works now.

gdalwarp -wo "SAMPLE_GRID=NO" -et 0 -rpc -to RPC_DEM=DEM_WSG84_ellipsoid.tif input.tif output.tif

I still get some small artifacts that disappear when I use -r bilinear but then I lose some contrast, so I am still trying to figure out how to improve this. But overall the results are almost exactly what I get when using dedicated GIS software.